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FOREWORD 

A  finite  difference  method  suitable  for  the  design  of  finned  bodies  in 
supersonic  flight  is  described.  Efficient  numerical  calculations  are  ac.hieved 
using  a  thin  fin  approximation  which  neglects  fin  thickness  but  retains  a 
description  of  the  fin  surface  slope.  The  resulting  algorithm  is  suitable 
for  treating  relatively  thin,  straight  fins,  with  sharp  edges  which  may  be 
deflected.  Methods  for  treating  the  fin  leading  and  trailing  edges  are  described 
which  are  dependent  on  the  Mach  number  of  the  flow  normal  to  the  leading  edge. 

The  leading  and  trailing  edge  analysis  varies  from  exact  to  empirical  as  the 
normal  component  varies  from  a  supersonic  to  subsonic.  A  procedure  for  modeling 
body  crossflow  separation  using  a  Kutta  condition  is  described  which  yields  a 
reasonable  leeside  vortex  structure.  Calculations  are  compared  to  experiment 
for  body-alone,  body-wing,  body-tail  and  body-wing- tail  configurations.  The 
computer  program,  developed  in  this  study,  is  described  in  a  separate  report. 
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SECTION  1 
INTRODUCTION 

A  practical  means  of  predicting  the  nonlinear,  inviscid,  supersonic  shock 
layer  on  missile  configurations  is  to  numerically  solve  the  steady, 
three-dimensional  inviscid  equations  using  an  efficient  finite  difference 
method.  Several  computer  programs  are  currently  available  for  calculating 
flow  fields  about  arbitrary  bodies  in  supersonic  flow.  However,  their  application 
to  practical  wing-body-tail  configurations  presents  some  serious  computational 
problems.  Existing  codes  treat  the  complete  fin-body  cross  section  as  a  single 
entity.  Thus  when  cylindrical  coordinates  described  in  Figure  1  are  used,  a  large 
number  of  ;*)  mesh  planes  are  needed  to  adequately  resolve  the  fin.  When 
several  fins  are  present  at  the  same  axial  station,  the  number  of  grid  points 
needed  becomes  prohibitively  large  for  practical  design  calculations.  The 
number  of  grid  points  can  be  substantially  reduced  by  mapping  the  fin  body 
cross-section  into  a  more  rounded  figure.  The  existing  methods  utilizing  this 
approach  are  based  on  conformal  mapping  techniques  developed  by  Moretti^’® 

(also  see  References  3  and  6).  However,  the  mappings  are  complicated  even  for 
the  case  of  a  single  smooth  fin  or  wing  and  often  tend  to  cluster  large  number 
of  mesh  points  near  wing  tips.  This  reduces  the  permissible  marching  step  and 
increases  computational  time. 


The  primary  focus  of  the  present  study  is  the  development  of  a  more  efficient 
numerical  technique  for  treating  finned  bodies.  To  achieve  this,  the  approach 
used  here  departs  from  the  basic  computational  strategy  used  in  References  1-7 
when  fin  surfaces  are  present.  Instead  of  considering  the  cross-sectional 
body-fin  geometry  as  a  single  entity,  the  present  method  considers  the  body 
alone  (i.e.,  the  body  with  all  fin  surfaces  removed)  and  the  fin  geometries 
separately.  The  computational  grid  is  generated  using  normalizing 
transformat ions^ ^ ^  applied  to  the  body  alone  co:.f iguration .  The  fin 
surfaces  are  allowed  to  extend  into  the  computational  region  and  can  be 
adequately  resolved  within  a  relatively  coarse  computational  grid.  In  order  to 
treat  the  complex  flow  in  the  immediate  vicinity  of  fin  leading  and  trailing 
edges,  appropriate  local  analyses  are  built  into  the  program  which  depend 
strongly  on  the  local  Mach  number  of  the  flow  component  normal  to  the  edge. 

These  local  analyses  can  range  from  locally  exact,  when  the  edge  is  sharp  and 
the  normal  velocity  component  is  sufficiently  supersonic,  to  hoc  or 
semi-empirical  in  other  situations.  It  is  possible  to  exercise  the  above 
computational  procedure  without  recourse  to  a  special  leading  edge  analysis. 
However,  such  a  procedure  is  not  as  robust  and  does  not  resolve  the  leading  edge 
region  as  accurately. 
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Within  this  framework,  various  approaches  for  numerically  treating  general 
fin  surface  shapes  are  possible.  One  approach  is  to  introduce  extra 
computational  points  to  represent  the  fin  surfaces  which  would  float  within  the 
basic  grid.  This  would  complicate  the  application  of  the  boundary  conditions  on 
the  fin  surfaces.  Another  approach,  is  to  subdivide  the  flow  domain  into 
several  sub-regions  each  containing  the  flow  between  adjacent  fin  surfaces. 
Relatively  simple  transformations  are  applied  separately  in  each  sub-region  to 
map  adjacent  fin  surfaces  onto  constant  computational  coordinate  planes. 
Relatively  coarse  meshes  could  be  used  in  each  sub-region  and  the  computations 
In  the  various  sub-regions  could  be  linked  in  a  manner  suggested  by  Hindman,  et 
al.^  Both  the  above  mentioned  approaches  are  in  principle  capable  of  handling 
general  fin  surface  geometries. 

To  simplify  the  development  for  the  present  study,  the  analysis  is 
restricted  to  relatively  thin  fins  with  sharp  edges  which  lie  approximately 
along  constant  $  planes.  A  thin  fin  approximation  is  employed  which  neglects 
the  fin  thickness  but  retains  the  actual  fin  surface  slopes.  For  an  important 
class  of  body-fin  configurations,  the  thin  fin  approximation  allows  the  direct 
use  of  the  basic  grid  generated  for  the  body  alone  shape  without  the 
introduction  of  floating  points  to  describe  the  fin  surface  or  additional 
mappings . 

A  computer  program  has  been  developed  using  the  thin  fin  algorithm  which  is 
capable  of  treating  configurations  with  a  large  number  of  lifting  surfaces.  The 
types  of  geometries  which  can  be  handled  are  restricted  as  follows: 

1 .  The  body  alone  surface,  b(2,<ti),  must  be  single  valued  in  'll.  Thi s 
precludes  a  direct  treatment  of  items  such  as  detached  inlets. 

2.  Fins  must  be  relatively  thin  and  lie  near  constant  ^  planes.  In 
practice  it  has  been  found  that  reasonable  agreement  can  be  obtained  between 
calculation  and  experiment  on  relatively  thick  or  deflected  fins.  By  moving  the 
coordinate  system  origin,  it  is  often  possible  to  position  fins  along  constant 

<:>  planes. 


3 .  Fin  edges  must  be  sharp  and  the  fin  edge  radial  location.  L(z),  can  be 
single  or  double  valued  in  z.  This  allows  swept  trailing  edge  configurations  to 
be  treated. 

4.  Fins  cannot  extend  through  the  bow  shock. 

In  addition  to  the  above  requirements,  the  flow  field  must  remain 
supersonic  throughout  the  entire  calculation.  This  precludes  bodies  with  blunt 
protuberances  which  feature  upstream  subsonic  flow. 

The  computational  method  described  in  this  report  is  an  extension  of  the 
method  developed  earlier  for  re-entry  bodies  with  cuts  and  flaps.  A  detailed 
description  of  this  method  is  contained  in  Reference  10  and  will  be  only  briefly 
outlined  here.  The  present  report  provides  a  description  of  the  analysis  for 
configurations  with  fins,  various  special  procedures  required  for  different 
configurations  and  flow  conditions ,  and  also  comparisons  of  the  numerical 
results  with  available  experimental  data.  A  description  of  the  computer 
program  and  user  guide  is  provided  in  Reference  11. 
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SECTION  2 

COMPUTATIONAL  ALGORITHM 


Euler's  equations  are  solved  numerically  in  a  region  bounded  by  the  missile 
body  (i.e.,  with  lifting  surfaces  excluded)  and  the  bow  shock.  At  interior 
points  these  equations  are  cast  in  weak  conservation  form  and  advanced  using  the 
MacCormack  predictor-corrector  method.  At  the  body,  fin  surfaces  and  bow  shock, 
special  sets  of  equations  are  solved  using  the  MacCormack  predictor-corrector 
method  with  appropriate  one-sided  differences.  These  sets  of  relations  combine 
the  admissible  characteristic  compatibility  relations  and  the  appropriate 
boundary  condition. 

2.1  COMPUTATIONAL  REGION 


The  flow  field  equations  are  initially  written  in  cylindrical  coordinates 
(r,'>,z),  shown  in  Figure  1,  and  then  transformed  into  (X,Y,Z)  coordinates  for 
the  purpose  of  numerical  integration.  The  current  computational  algorithm 
admits  transformations  of  the  form: 

Z  -  2,  X  •  X(r,«  ,2),  Y  =  Y(4>  ,2)  (1) 

This  transformation  maps  each  crossflow  plane  into  the  square  computational 
domain  depicted  in  Figure  2.  The  missile  body  (i.e.,  excluding  the  lifting 
surfaces)  and  shock  are  located  at  X-0  and  X“1  respectively  while  the  Y=0  and 
Y“1  planes  correspond  to  ♦  0  and  ♦  =•  ♦*.  Fin  surfaces  are  allowed  to 

extend  Into  the  computational  domain  along  constant  Y  planes.  The  thickness  of 
fins  is  neglected  via  the  thin  fin  assumption  which  is  described  later  in  this 
section.  Two  sets  of  variable  values  are  carried  at  grid  points  describing  a 
fin,  one  for  the  upper  surface  and  the  other  for  the  lower  surface.  The 
required  mapping  is  thus  only  a  function  of  the  body  alone  plus  shock  geometry 
and  it  is  expressed  as  a  composite  of  two  transformations.  The  first  is  the 
usual  normalizing  transformation  given  by 


X  =  [r-bC^  ,z)  J/(c(<i  ,z)-b(i(i  ,2)  ] 
y  2  »  2 


(2) 


The  second  mapping  is  primarily  used  to  cluster  computational  points  in  r,i>,z 
space  while  retaining  a  uniform  grid  in  the  computational  region.  This  mapping 
IS  conveniently  expressed  in  inverted  form  as 


X  -  f(X,Y,Z),  y  =«  g(Y,Z),  2  =  Z 


(3) 
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where  f(0,Y,Z)  =  0,  f(l,Y,Z)  »  1;  g(0,Z)  =  0,  g(l,Z)  -  1.  Apart  from  these 
conditions,  the  mapping  functions  f  and  g  are  arbitrary  provided  that  the 
functions  and  their  derivatives  through  the  second  order  are  smoothly  defined 
and  fx  and  gy  are  positive.  In  the  present  code  the  functions  f  and  g  can 
be  analytically  or  numerically  defined  as  is  discussed  in  Reference  11.  In  the 
latter  case  the  required  derivatives  are  defined  using  second  order  finite 
difference  methods.  If  mesh  clustering  is  not  required,  f  =  X  and  g  =  Y. 

From  (2)  and  (3): 

^2  “  “8z/8Y»  “  l/('t'*8Y). 

Xr  -  l/[fx(c-b)], 

(A) 

Xz  -  -(f2+Y2fY)/fx^■Xr[(f-l)b2-fc2]. 

X0  =  -Y<,fY/fx+Xr[(f-l)b,t>-fc,),  ]. 

When  there  are  no  fins  present  each  computational  plane,  Z  =  constant,  is 
covered  by  a  uniform  mesh  defined  by 

{(Xn,Yn,)  :  Xf,  =  (n-l)4X,  Yn,  =  (m-l)AY,  n  =  1,2. ..,N;  m  *  1,2. ..M) 

where  AX  »  1/(N-1)  and  AY  =  1/(M-1)  for  the  symmetric  case  and  AY  =  1/M 
for  the  asymmetric  case  (N,M  are  positive  Integers).  When  fins  are  present, 
this  computational  mesh  is  modified  by  applying  planar  cuts  along  the  fin 
surfaces . 

The  algorithm  for  advancing  the  unknown  flow  field  quantities  from  Z  = 
to  z'^'*'*'  =  Z^iZ  depends  on  the  location  of  the  individual  mesh  point  in  the  shock 
layer.  These  are  divided  into  the  following  four  types:  interior,  body  surface, 
shock,  and  fin  surface  points.  The  procedure  for  treating  the  first  three  types 
of  points  is  outlined  in  the  next  section  with  more  details  provided  in 
Reference  10.  This  is  followed  by  a  discussion  of  the  fin  surface  treatment.  The 
derivation  of  the  fin  surface  compatibility  relations  is  provided  in  Appendix  A. 

2.2  INTERIOR,  BODY,  AND  SHOCK  POINTS 


The  flow  field  quantities  at  interior  points  are  advanced  using  Euler’s 
equations  in  weak  conservation  form: 

3  (rU)  +  3  (rF)  +90=  'e  (5) 

3  z  3  r  34 


where  U,F,G,  and  E  are  column  vectors  defined  in  transpose  form  by 

=  (Pw,  p-tpw2,  owu,  pwv) 

F^  =  (pu,  pwu,  p+Pu^,  Pvu) 

C.t  =  (pv,  p  vw,  pvu,  p-tpv^) 

Et  =  (0,  0,  p+Pv2,  -puv). 


12 


NSWC  TR  81-457 


The  energy  equation  for  a  steady  inviscid  flow  with  an  isoenergetic  free 
stream  reduces  to  the  following  algebraic  relation: 

h  +  (u2+v2+w2)  =  h  +i(u2+v2+w2)  *  const*  =  H  (7) 

7  00  2  00  00  00  O 

The  gas  is  assumed  to  be  perfect  which  allows  the  system  to  be  closed  by 
applying  the  definition: 

h  -  VP  (8a) 

(Y  -1  )P 

Additional  thermodynamic  relations  which  are  used  throughout  this  report  are: 

a^  =  Y  p/p  (8b) 

s  ”  In  p  -  Yin  p  (8c) 

The  system  (5)  when  transformed  by  (1)  can  be  written  as 

^  +  y;  +  ^  .  E  (9) 

az  ax  ay 


where 


U  =•  rJ-lu,  F  =  rJ-1  [XzU+X^?+(X^/r)5')  , 

G  =  rJ-l[YzlJ+(Y^/r)G],  E  =  (J-l)E,  J  =  (X^Y^) 


The  above  equations  are  integrated  using  MacCormack  differencing: 


U* 


uk+1  =  1 
n,m  2 


^tt+I,m  ~ 

^  rrt-I-l,m 

1  AZ  -| 

r 

^n,m+J  ~  ^n,nH-J' 

AX  J 

AY 

n,nH-J-l  Vz  +  E^  „AZ  + 
- 1  -n,m 


^^n,m''’  ^n,m  ”  ^ 


(10a) 


F  -F*\  /g*  -G*  \  *1 

tr4-l-I,m  n-I>c^  AZ  -  ^  n.m+l-J  n,m-J j  az  + 

(10b) 


AX 


.  k  +  1 

where  U*  is  the  predictor  value  and  F*  =  F(U*,Z  ,  X  , -Y  ). 

n,!n  m  m 

Here  I,  J  =  1  produces  a  forward  difference  in  the  predictor  step  and  a  backward 
difference  in  the  corrector  while  I,  J  =•  0  yields  a  backward  difference  in  the 
predictor  step  and  forward  difference  in  the  corrector  step*  At  the  end  of  both 
the  predictor  and  corrector  steps  the  physical  variables  p,P  ,u,v,w  are  decoded 
from  U  and  (7)  using  the  relations: 


w 


[Y  +  ^1  -  ] 

Ui  ( Y  +  1) 


(y2  -  1)  [H 

o 


U 


i)2 

U2 


I] 


(11) 
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u  -  U3/U1,  V  =  U4/U1 
P  -  (L’2  -  Uj^w)(J/r);  D  =  (Ui/w)  (J/r) 

The  step  size,  AZ,  is  chosen  to  satisfy  the  CFL  stability  conditions  for 
the  MacCormack  scheme.  This  consists  of  ensuring  that  the  domain  of  dependence 
of  the  partial  differential  equations  Is  contained  within  the  domain  of 
dependence  of  the  finite  difference  equation  at  all  points.  The  derivation  of 
the  CFL  condition  for  the  MacCormack  scheme  Is  provided  In  Reference  10.  The 
resulting  relations  are: 


where : 


4Z  =«  c  AX  min  {  } 

U 

M  =•  max(ui,y  2»^*  3) 

U  =  jwA-a^X^j  +  a  ^w2-a2)(x2  +  X|/r2)  +  (A-wXp2 

^2  *  1^1  a  y(w2+v2-a2)(Y^r2)  ] 

=  [wA-a2x2  -  ^(wB-a2Y2)| 


(12) 


%  (  AX/AY,  if 

6  ”  ) 

(  -  AX/AY,  If 


I  -  J 
1  ¥  J 


A 


X  u  +  V  +  X  w; 
r  r  z 


»  +  Y  w 

r  z 


5  =  safety  factor  .9 


The  minimum  iZ  taken  over  all  the  computational  points  is  used  as  the  step  size 

On  the  body  boundary  (X  =  0)  the  component  of  velocity  normal  to  the  body 
surface  must  vanish;  l.e.. 


u  -  bjW  -  (b^/b)v  =  0.  (13) 

This  condition  is  supplemented  with  certain  characteristic  compatibility 
relations  associated  with  the  system  (9).  This  approach  was  originally 
suggested  by  Kentzerl2  unsteady  flow  problems  and  has  been  adapted,  in 
different  form,  to  steady  supersonic  flow  In  References  4  and  5.  The  method 
used  In  this  work  is  given  in  detail  in  Reference  10  and  is  only  briefly 
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reviewed.  It  is  found  that  there  are  three  independent  characteristic  relations 
which  are  admissible  on  X  >  0.  These  can  be  written  as  a  system  of  quasl-llnear 
first  order  partial  differential  equations  on  X”0  for  advancing  P  “  ln(p), 

V2  =  u  ^  +  V  and  s.  The  resulting  relations  are: 

IZ  ”  fx  X  1  {  pw[  X  -  (a  w  +  a  v)]  +  ^1  i 

3Z  L  r  -MX  ^  ^  3X  7  4  J  p  (14a 

A 

where  p  ■  £21  V  +  wX  e  -t-  C  *  (14b 

b  2  -t-  1  3Y 


a2  (/3-b  ) 

^  - - L-£_  .  0 

+  22  1 
w  -a 


w_V2  -  [1  +  cj_)  2| 


‘  -  hz<^Yz/Y^ 

/  9  Z 


s(K/b)  ,  b  b.  , 

•4  --si—  -  <»>»♦  - 


?  =  (5l,C2t5  3,54) 

C  *  wA  (2-V^  /p)  ,C  =«b-X  +  w^  X  /p 

1+  1  2z  +  l-(- 

C  =•  wu<  X  /p -1  ,  5  “  wvK  X  /p  -H  b  /b 

3  1  4  1  -t-  1. 

“  (1^)  “  1/ (2JL)  (<  “  -  p/h  for  a  perfect  gas) 

1  3h  p  3p  p  1 

^85  "  “  8Zy/8Y  .  '2  =  GJ/r 


=  1  +  +  b2 


a  u  -t-  "/p  w 
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where 

A 

-fr-|G  -  pvwb^/b  -|(Tg5Y^  +  b^/b) 


(15b) 


7  -  (’)i,’?2.’»3.’J4) 

ni  =  V2  ,  02  •  0  ,  03  «  -  b|j,/b  ,  04  “  -1 


3_s  =  -  B  ii  (16) 

az  w  3Y 

Equations  (14)  to  (16)  are  advanced  using  a  predictor-corrector  method  of 
the  form  (10)  but  with  the  X  derivatives  replaced  by  forward  differences  in  both 

the  predictor  and  corrector  steps.  This  differencing  scheme  which  is  first 

order  accurate  can  be  made  second  order  accurate  using  the  procedure  given  in 
Reference  10.  At  the  end  of  both  the  predictor  and  corrector  steps,  wall  values 
of  p,p,u,v,w  are  determined  from  P,V2,s, (7) , (8c)  and  (13)  which  can  be  manipulated 
to  give  the  relations: 

p  =  exp(P);  q2  =  2(Ho  -  h) 
p  =  (p/exp(s))l/Y 


w 


Vfi+<^/b)2]q2-v2  /  vy 


(17) 


V  =  [V2  -  (b^/b)b^w]/[l  +  (b^/b)2] 
u  =  bzW  +  (i^ji  /b)v  . 

A  A 

Alternative  expressions  for  p  and  v  of  (14b)  and  (15b)  which  often  give 
improved  results  are: 

p  -  ^  +  [b.Y^  b*Y^  +  (Bw  -  Y^)  X+]|i  (14c) 

b  a 


-  oB  (a5W  +  a3v) 

+  PwA  {  (Tf  -  T„  )  w  +  Z  Tf  +  WfY.  ^  (v/w)  +  b  1} 
+  ^6  ®6  b  ‘7  b  ^  aY  ^ 
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where 


and 


3  b 

»  - _ I 

5  3  Y 


/b) 

a,  -  — I—  -  [b^^/b  -  (b^/b)2]/Y^ 

^  0  1 


f  +  Y  f  (b  -  c  ) 

T.  -  T  +  ZX  2  YX  -  z  2 

f6  86  - - -  - - n — 

fx  (c  -  b) 


fx  (c  -  b) 


3  V 

V  "PBCa^u  - 

3  3  Y 


2)  -i  y1£_  -  pvw  b 
b  $  3Y  ~5“  2 


(15c) 


These  equations  are  algebraically  identical  to  (14b)  and  (15b).  They 
produce  slightly  different  numerical  results  because  the  X  and  Y  differences 
Involve  different  quantities. 

Many  body  configurations  of  Interest  have  sharp  corners  or  edges  such  as 
those  found  on  blconlcs  and  other  segmented  shapes.  If  the  upstream  body 
surface  velocity  normal  to  this  edge  Is  supersonic,  either  a  shock  wave  or  an 
expansion  fan  will  be  attached  to  this  edge  producing  a  discontinuity  In  surface 
flow  variables.  To  handle  this  situation,  an  oblique  shock  or  Prandtl-Meyer 
expansion  Is  applied  at  the  edge  as  Is  described  In  References  8  and  10.  In  the 
interior  these  discontinuities  are  captured  using  the  dissipative  and 
conservation  properties  of  the  Interior  point  scheme. 


At  the  bow  shock,  flow  field  properties,  as  well  as  the  shock  slope,  are 
unknown.  The  correct  boundary  conditions  are  provided  by  the  Rankine-Hugonlot 
conditions  which  relate  the  free  stream  properties,  the  shock  slopes,  and  the 
properties  behind  the  shock.  An  analysis  of  the  characteristics  associated  with 
the  system  (9)  indicates  there  Is  one  admissible  characteristic  rel.-ion  on  X  =  1 
(see  Reference  10) .  This  expression,  when  combined  with  the  Ranki;  e-Hugoniot 
relations,  results  in  the  following  system  of  equations  which  are  usel  to 
advance  c,  C2  and  c,^  : 


3  c  « 

c  -  (Y  /Y  )  c 

3Z 

z 

z  <p  Ip 

3  c . 

3c 

3  c 

♦ 

a  Y  _ Z  -  Y  _ z. 

3  Z 

♦  3Y 

z  3Y 

C.  3  c 

3  c 

c 

1 

L 

L 

-  Y  _ 5.7  -  .»  (r 

s  c  3  Y 

2  3  Y 

c 

Y  c  /Y  'I  ]} 

z  ^  ^ 


(18) 
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:  (Vg  w«  -  Cj  Vji  )  Ai  -  (p  -  p»)[l  +  (c,>/c)  ]J/'' 


:[Vs  V,  -  (c^/c)Vj^]  Ai  +  cz(c^/c)(p  -  f^))vg 


[  +  P  s  ^  ~  Cz  Vjj  )  ]  Aq  +  Aq  Poo  (Vf 


^Vn  -  Vn)[a^  +  Vn  +<i  (a^P)V„(V„  -V.)] 


±  \  -  a^)[l  +  (c  /c)^]  +  [u  -  (c  /c)v]^ 

a  T  <>  ♦ 

(i£.)  “  (•«  ®  ~  £.  for  a  perfect  gas) 

3h  p  3p  p  1  h 


A  A  A 

c*  (i£  +  12  +  E)  -  (A  -  V  pw){v  +  (c  /c)u  ]Y  /(Y  v  ) 
3X3Y  Is  *  ^®zi>s 


3X  3Y  ~3X  -  3X 

FJ/r  ;  E  =(F  -  E)-  f _ £  + _ 2]  U  - _ E  F  -  [ _ i 

- r —  3X  3Y  3X  3X 


?2»  5  3»  ?4) 

2 

[2  -  (ki/p  )V  ]X_,  ?  2  =  [(u  ~  (c,^/c)v)  -  X_]/w  +  Ki 

f  1  u/p  -1,  C4=c,j,/c+X_ici  v/p 

2  2  2 
-a  {  AqW  +  [u  -  (c^/c)vj}/(w  -  a  ) 


-  C  _ 

(u  -  VC  /c  -  wc  )/v  ;  v4  »  1  +(_2.)  + 
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In  the  above,  the  quantity  C**  (LL  +  L2.  +  E)  denotes  the  inner  product 

3X  3  Y 

of  these  vectors  and  Vjj  is  the  free  stream  velocity  component  normal  to  the 

00 

shock.  Since  Cj  >  0,  the  equations  (15)  -  (17)  can  be  used  to  advance  c, 
c^,  and  c^  in  a  predictor-corrector  method  of  the  form  (10)  but  with  the  X 
derivatives  replaced  by  backward  differences  in  both  the  predictor  and  corrector 
steps.  With  the  local  shock  angles  thus  determined,  the  flow  variables  at  X  =  1 
follow  from  the  Rankine-Hugoniot  conditions  and  tne  known  free-stream 
conditions.  The  resulting  relations  are: 

p  •  _i_  fp  (1-^)  +  2  p  v2  I 
Y+1  ®  »  Ha  J 

(19) 

2  2  2 

P  -  P  Vn  /(Po,  Vn  +  pto  -  p) 


0  \ 

_!L  J,W®W  -(u-u)c,V‘»V  -(u-u  )(c  /c) 

P  /  m  o>Z  »  a>^ 

In  the  above,  Ho,^to,Ho  are  the  free  stream  velocity  components  given 
by  Voo  cosfl  cosa ,  V*  (cosfl  sin  a  sin^  -  sinB  cosi^),  and 
-V^  (cosB  sina  cos<t>  +  sinB  sin(ti),  respectively. 

2.3  FIN  SURFACES 

2.3.1  THE  THIN  FIN  APPROXIMATION  AND  IMPLEMENTATION.  The  thin  fin 
approximation  is  applicable  to  fins  with  surfaces  that  lie  close  to  a  constant 
♦  plane,  say  ♦  ®  which  is  defined  as  the  fin  plane.  The  fin  geometry 
is  assumed  to  be  represented  by  two  surfaces,  the  upper  and  lower  surfaces,  each 
described  independently  by  relations  of  the  form 

$  =  <^f  +  a(r,z)  (20) 

In  the  cross-section  Z  •  constant,  the  actual  fin  surfaces  will  lie  within  the 
computational  mesh  as  shown  in  Figure  2.  The  thin  fin  approximation  assumes 
that  o  is  small  and  thus  places  the  fin  surfaces  along  the  fin  plane 
corresponding  to  Y  -  Yf  in  each  Z  =  constant  plane.  Although  the  fin  is 
approximated  by  a  zero  thickness  plane  lying  on  ^  “  ♦f,  the  surface  slopes 
are  described  to  0(|aj).  The  fin  surface  is  prescribed  by  specifying 
0 (r,z)  ,v(r,z) ,  L(z)  and  the  first  derivatives  of  these  quantities.  Here 
9  and  v  are  the  angles  between  the  fin  surface  tangency  plane  and  the  fin 
planes  in  the  r  and  z  directions  respectively,  and  the  quantity  L  is  the 
radial  location  of  the  fin  edge.  In  terms  of  9 ,  v  the  derivatives  of  a 
correct  to  0(|cj  )  are  given  by: 

ro  “  tan  9  ,  ra  ^  =  tan  v  , 

ra  J.J.  »  sec^  9(9j.-i;j.)-ar  (21) 
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^  1  i 

r«2z  "  sec  vv^  -  \  tan  v,  r'7  “  sec  9  (6  ^  -  a 

Within  the  restriction  that  |  o|  be  "small",  the  thin  fin  approximation  can  be 

applied  to  arbitrary  fin  geometries  including  surfaces  with  discontinuous  slopes 
and  fins  with  "small"  deflections,  camber,  and  variations  in  dihedral. 

The  numerical  algorithm  for  treating  fins  by  the  thin  fin  approximation 
requires  that  the  computational  mesh  be  chosen  so  that  each  fin  plane  is 
coincident  with  a  computational  mesh  plane,  Y  =  Yf.  Two  sets  of  computational 
points  are  carried  on  the  Y  =  Yf  plane  to  describe  the  flow  properties  on  the 
upper  and  lower  surfaces  as  is  illustrated  in  Figure  2.  The  upper  surface  of  a 

fin  features  an  outward  normal  that  has  a  positive  ^  component,  while  the 

lower  surface  has  a  negative  ^  component  associated  with  its  outward 
normal.  As  the  calculation  is  marched  down  the  length  of  the  body,  fin  surfaces 
are  encountered  on  Y  **  Yf.  Thus  a  point  at  some  X  may  at  one  axial  location 
be  an  interior  flow  field  point  and  in  the  next  axial  step  move  onto  the  fin.  Here 
the  Interior  point  is  split  into  two  points  corresponding  to  the  upper  and  lower 
fin  surfaces.  The  fin  points  thus  created  are  referred  to  as  leading  edge 
points .  For  a  fixed  X,  a  pair  of  points  which  are  on  the  fin  at  one  axial  step 
can  in  the  next  step  move  off  the  fin  and  become  a  single  interior  flow  field 
point.  Such  a  point  will  be  referred  to  as  a  trailing  edge  point.  The  flow 
variables  at  leading  and  trailing  edge  points  are  determined  from  an  appropriate 
local  analysis  which  is  described  in  the  following  subsections.  The  adjustment 
for  the  presence  of  a  leading  or  trailing  edge  is  made  immediately  after  the 
completion  of  the  step  in  which  the  edge  is  encountered.  The  values  of  the  flow 
variables  prior  to  the  adjustment  are  termed  upstream  while  the  adjusted  values 
are  termed  downstream.  Note  that  the  locations  of  leading  and  trailing  edge 
points  are  within  one  42  of  the  physical  edges  of  the  wing. 

2.3.2  FIN  SURFACE  BOUNDARY  CONDITIONS.  On  a  fin  surface,  the  velocity 
component  normal  to  the  surface  must  vanish;  i.e., 

v/r~<y^w-<j^u“0.  (22) 

The  numerical  methods  used  to  advance  the  fin  surface  points  are  based  on  the 
appropriate  charactistlc  compatibility  relations  associated  with  (9)  which  are 
derived  in  Appendix  A.  Both  the  upper  and  lower  fin  surfaces,  although 
considered  separately,  are  treated  using  the  same  techniques.  For  fin  surface 
points  not  on  the  fin  body  junction  (X=0),  the  three  compatibility  relations 
listed  below  are  used  to  advance  s,  V3  =  u  +  roj-v  and  P  =  in  p  along  the 
fin  surface: 

Li  =  -  i(X  w  +  Xu  +  X.v/r)Li 


=  v(a2  +  a3)  +  V/pw 


3(X^/r). 


mere  V  =  (-V  ,0,l,ro^)  K  F  )  +  p[ _ r  +  ^o 

3  3X  ax  "  ax 


rz  -  («r  +  tnrr)ao.  «0  ~  /  (Xj.  +  c^X^) 
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=  Xy^  +  Dw(a2  ^  +  Oj.  ^  -  1  8v)  ] 

9Z  3Y  S  3Y  3Y  r  3Y 

pw[u(a]^  -  X)  +  wa2  +  vV2/(rw)]/r 

-rr  3'F  +  (pw2x  +  pri9)^^2  +  (pwuX  +  pri9)^^r 

3X  3X  3X 

+  (pwvX  +  pn/^)^^X(j)/r) 

3X 

where  3  =  +  [w2(i  +  ai  +  o?)  -  1  - 

■  II?2  72 

A  =  a2(3  -  Cz)  ,  Yr  =  Y^X^/CX^  +  a^X^) 
w2-a^ 

2 

*2  =  ”22  -  “o(°2  +  Wrg),  ni  =  Xw(2  -  <  |q|  ), 

2 

n2  “  (w  <  “  1)X  +  <*21  »i3  =  wtjXk  +  o  r> 

n4  •>  wvXie  -  1/r 

0K  "  (3p/3h)  =-  £  for  a  perfect  gas 

P  h 

The  upper  and  lower  signs  in  B  are  used  on  the  upper  and  lower  fin  surfaces, 

respectively.  The  quantities  and  in  (23)  (24),  and  (25)  represent  the 

3  Z  3  X 

partial  derivatives  of  quantities  defined  on  the  fin  surface  and  thus  are 
functions  of  Z  and  X  only.  Equations  (23),  (24),  (25)  follow  from  the 
compatibility  equations  derived  in  Appendix  A  without  approximation.  The  thin  fin 
assumption  is  applied  by  evaluating  flow  field  variables  and  transformation 
quantities  (Xy ,X2,J^  ,etc. , )  at  the  fin  plane  rather  than  at  the  actual  fin 
surface.  Predictor-corrector  differencing  of  the  form  of  (10)  is  used  to  march 
the  solution  in  Z.  However,  the  Y  derivatives  are  replaced  by  forward  differences 
on  the  upper  fin  surface  and  the  backward  ones  on  the  lower  fin  surface  in  both 
the  predictor  and  corrector  steps.  At  the  end  of  each  computational  step,  the 

quantities  p,p,u,v,w  are  determined  by  simultaneously  solving  (7),(8c)(22)  and 
the  definition  of  V3.  This  yields: 

p  »  exp  (P)  ^  (26a) 

P  =  [p/exp(s)]  V  (26b) 


« 
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/2(Ho  -  h)(l  + 


'22^  2,2  ,  , 
r  (Jj.  +  r  ’2  +  1 


(26c) 


(•to  \V.  +  ro  w 
V  s'"  r-'  3 _ 2 


(26d) 


1  +  r"  0: 


u  =  V3  -  10  j-V 


(26e) 


The  fin  body  junction  is  assumed  to  be  a  sharp  corner.  At  this  corner  both 
(13)  and  (22)  are  satisfied  and  thus  flow  is  directed  along  the  corner.  This 
implies  that  entropy  is  constant  along  the  corner  except  at  compressive 
discontinuities  in  body  or  fin  slope  and  at  leading  and  trailing  edges.  Since 
(7)  also  holds,  only  one  additional  relation  is  needed  to  completely  determine 
the  flow  variables  along  the  junction.  This  is  given  by  a  characteristic 
compatibility  condition.  However,  an  ambiguity  arises  in  the  choice  of  an 
appropriate  characteristic  condition.  Two  possible  equations  are  derived  in 
Appendix  A: 


i£  ”  X_A  1  [Le  +  iii(  b  i”  +  JL  “  Lii)  ] 

az.  -^ax  132  b3X  3X 


+  £i£_  {  bYrAi  (a_  Le  +  a_  iw  -  1,  3V) 
bfli  ^3Y  3Y  b3Y 

+  w  (a^  +  02  A2)  +u(a2A2  ~A2)  +^[05 
-  (V2  -  A2V3/b)/w]} 


(27a) 


where 


Ai  *  a  (a2T  -  b2R  +  fliR)AJi,  A 2  =  (-«2 


2  2  2  2  2  1/2 
+  a  Qi  =  {  ~  n2)/[a  (w  "  a  )R]} 

2  2  2  2  2  2 
R-02+0r  +  l'^b,S*bz  +  l  +  h),/b  ,T=a2bz 

2  2  2  2  2  2 
-Oi.-b^/b,ni“(w  -  a  )R  +  a  a  ^  2  ~ 

2  2  2  2  2  2 
-a  )T  +  a  b^z»  ^  3  *  (w  -  a  )S  +  a  b^, 

'»4  =•  fc(fcz2  'ibaji),  a  3  =  (cr  fbz  +  o  2)/(l  -  c  ^b^  )  , 


‘5  “  -  h^bg/b  +  a  3(bb^,>  -  b^,  )/b  . 
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and 

LE  •  Y  A  [LE  +  L2(o  iZ  +  a  LE-A  iZ)  I 
3Z  C33Y  fl2  z3Y  r3Y  b3Y 


b 

+  pw  (  bXA  (b  A5i-LE  +  _!.  iz) 
bfl2  r323X  3X  b3X 


+  w(o^^  +  02)  +  u(a2  -  A3)  +  v[A4a5 

-  (A4V2  -  V3/b)/w]}  (27b) 


2 

where  A3  “  a  [  b^T  -  a  gS  +  Sfl2]/Q3,  A 4  *  ^"^2 

2  1/2 
+  a  b2fl2)'^3>  ®2  "  i  fli(R/S)  .  The  upper  and  lower 

signs  are  used  for  the  upper  and  lower  surface  junctions,  respectively.  In 

(27),  3_  and  3 _  represent  the  derivatives  of  quantities  defined  on  the  fin 

3X  3Z 

surface  as  a  function  of  X  and  Z  only.  Equation  (27a)  is  based  on  the 
bicharacteristic  lying  along  the  wing  surface  while  (27b)  corresponds  to  the 
blcharacteristlc  lying  along  the  body  surface.  Both  equations  appear  to  produce 
similar  results  except  when  large  pressure  gradients  occur  in  the  vicinity  of 
the  junction  (l.e.,  close  to  a  leading  edge  or  surface  discontinuities). 
Computational  experience  seems  to  Indicate  that  the  more  robust  relation 
corresponds  to  the  equation  lying  al  'ng  the  direction  with  the  smallest  pressure 
gradients.  Both  of  the  above  equa'xons  follow  from  the  compatibility  equations 
derived  without  approximation.  The  thin  fin  assumption  is  introduced  by 
evaluating  the  flow  field  variables  and  transformation  quantities  at  the  fin 
plane  rather  than  at  the  actual  fin  surface. 

At  the  completion  of  the  predictor-corrector  sequence,  the  pressure  along 
the  corner  has  been  determined  using  equations  (27a)  or  (27b)  and  the  entropy  is 
unchanged  fom  the  value  at  the  previous  step.  The  quantity  p  is  obtained  from 
(8c),  the  enthalpy  from  (8a)  and  the  magnitude  of  the  velocity  from  (7).  The 
velocity  components  u,v,w  are  resolved  by  requiring  that  the  velocity  vector  be 
coincident  with  the  corner  direction  which  is  given  by  iff  X  itb,  where  iff  and 
n^  are  the  normal  vectors  to  the  fin  and  body  surface  respectively,  as  shown 
in  Figure  3.  Hence, 


( b  +  a  b  ) 

|q|; 


b(a  b  +  P  ) 
V  =  r  z  z  I 


q! 


w  = _ z  » 


where 


I  "cl  '  “  Ind  ’  '  l^cl 

ifc  “  (bz  +  P  zb^  )er  +  b(p  ^bz  +  p  z)e^  +  (l"®r  b^  )ez 


(28) 
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2.4  FIM  EDGES  AND  SURFACE  SLOPE  DISCONTINUITIES 


2.4.1  LEADING  EDGES.  A  leading  edge  point  occurs  when  an  Interior  point 
moves  onto  the  fin  surface.  This  point  must  then  be  split  into  two  points  in 
order  to  allow  both  the  upper  and  lower  fin  surfaces  to  be  described.  Several 
different  strategies  are  available  for  treating  leading  edge  points.  The 
simplest  approach  is  to  switch  from  the  interior  point  advancement  scheme  (9)  to 
the  fin  surface  scheme  (23),  (24),  (25)  at  the  leading  edge.  This  is  referred 
to  as  option  0.  To  accomplish  this,  the  computation  is  carried  out  through  the 
predictor  step  using  the  interior  point  advancement  scheme.  At  the  end  of  the 
predictor  step  the  z  value  is  advanced,  the  fin  geometry  is  updated,  and  the 
leading  edge  is  detected.  The  corrector  step  at  the  leading  edge  point  is 
completed  using  the  fin  surface  compatibility  equations.  Thus  for  a  leading 
edge  occurlng  in  step  k  at  (Xj,,?^,): 


pk 

=  (p^“i  + 

■  p* 

.  AZ)/2 

n  ,!n 

n  ,m 

n,m  dZ 

1  n. 

dV- 

Vk 

=  (V^-l 

+  V*  + 

3 

A2)/2 

(29) 

3 

3 

3 

dZ 

n  »ni 

n,in 

n  ,m 

n,m 

sk  s 

(sk“l  + 

s*  +  Al 

^Z)I2 

n  ,m 

n  ,m 

n  ,m  dZ 

n,m 

Here  P*,V3  and  s*  are  the  values  at  the  end  of  the  predictor  sequence. 

To  evaluate  the  above  requires  that  P*,V3,s*  and  s'^"!  be 

constructed.  The  quantities  P  ,P,s  and  s*  follow  directly  from  calculated  p*.p 
y*  and  p  values.  To  form  Vj  =  u  +  (rCj.)v  and  V3*=  u*  +  (rc^)v*  the  value  of 
(rOj.)  at  the  leading  edge  point  is  used.  Leading  edges  at  the  body  fin 
junction  are  treated  in  a  similar  fashion,  but  in  this  case  only  the  pressure  is 
advanced . 


The  leading  edge  option  0  represents  a  formal  descretlzation  of  the  various 
applicable  equations  without  recourse  to  additional  modeling  at  the  leading 
edge.  In  principle,  for  this  option  to  be  successful,  sufficiently  fine  grid 
and  step  size  must  be  used  to  allow  the  solution  to  "capture"  the  effects 
associated  with  the  presence  of  the  fin  edge.  In  practice,  this  option  gives 
poor  results  or  fails  altogether  when  the  flow  features  shock  or  expansion  waves 
which  produce  large  jumps  in  the  vicinity  of  the  edge. 

An  alternative  approach  is  to  apply  locally  an  analysis  which  models  the 
flow  very  near  to  the  leading  edge.  This  is  designated  as  option  1.  The 
justification  for  this  option  is  that  in  most  calculations,  for  reasons  of 
computational  efficiency,  the  mesh  spacing  in  the  vicinity  of  the  edge  is  not 
sufficiently  fine  for  option  0  to  yield  satisfactory  results.  The  computational 
algorithm  proceeds  by  completing  the  step  in  which  a  leading  edge  is  encountered 
without  taking  the  fin  surface  into  account.  The  resulting  flox-?  properties  are 
then  taken  as  the  conditions  immediately  upstream  of  the  leading  edge.  An 
appropriate  local  analysis  is  then  used  to  determine  the  flow  quantities 
immediately  downstream  of  the  edge  for  both  the  upper  and  lower  fin  surfaces. 

The  downstream  flow  quantities  are  then  assigned  to  the  leading  edge  points  on 
the  appropriate  fin  surface.  The  fin  surface  variables ,P ,V3 , s  are  then 
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advanced  from  the  downstream  values  using  the  procedures  described  In  Section 

2.3. 

The  velocity  vectors  used  to  characterize  the  local  flow  at  the  leading 
edge  are  shown  in  Figure  3  where  nf  Is  the  fin  surface  normal  and  n_  Is  the 
normal  to  the  plane  defined  by  the  upstream  velocity  vector  q  and  leading  edge 
tangent  vector  t.  These  vectors  are  defined  by: 

T  ez 

nf  =  -  r  0  ^  -  r  o  z?2  ( ^0) 

'n_  =  vcj.  +  -  u)e^-  Lzve^ 

Depending  on  the  orientations  of  nf and  q,  the  local  analysis  introduces  an 
appropriate  expansion  or  compression  turn  of  q  from  the  plane  normal  to  n_  into 
the  plane  normal  to  nf.  The  specific  turning  procedure,  which  strongly 
depends  on  (the  Mach  number  of  the  upstream  velocity  component  normal  to 
the  edge)  is  described  in  Section  2.4.4.  The  local  analysis  is  exact  If  the 
turn  on  both  sides  of  the  wing  can  be  accomplished  either  by  a  Prandtl-Meyer 
expansion  or  an  oblique  shock  attached  to  the  leading  edge.  In  other  cases,  the 
concept  of  the  exact  local  analysis  breaks  down  and  approximate  or  heuristic 
local  analysis  procedures  are  used  as  outlined  in  Section  2.4.4. 

In  cases  when  Option  1  is  used  and  an  attached  shock  does  not  exist  at  the 
leading  edge,  it  has  been  found  necessary  to  specify  the  downstream  streamline 
direction.  Figure  4  illustrates  the  influence  of  changes  in  streamline  direction 
on  the  calculated  surface  pressures.  As  the  streamline  direction  is  turned 
outward,  the  pressure  gradient  downstream  of  the  leading  edge  becomes  increasingly 
negative.  In  the  current  method  the  streamline  direction  is  set  such  that 
tan  (u/w)  =  .09911. 

At  the  body  fin  junction  option  1  is  implemented  by  rotating  the  leading  edge 
flow  velocity  vector  within  the  tangency  plane  until  it  is  in  the  corner  direction, 
n^,  which  is  given  by  (28).  The  required  turning  angle  is: 

«  =  +  sin~l[(q-nf)/(  |q  l|nf  I)  ]  (31) 

The  -  and  +  sign  apply  to  the  upper  and  lower  surface  respectively.  An  expansion 
occurs  if  S  <  0  and  a  compression  if  6  >  0.  The  turning  procedures  described  in 
Section  2.4.4  are  used  to  determine  the  downstream  values  of  p,  P ,  and  q.  It 
Is  possible  for  a  body  slope  discontinuity  to  occur  at  the  same  axial  station  as 
a  body-fin  junction  leading  edge.  In  this  case  the  body  discontinuity  Is 
treated  first  with  the  presence  of  the  fin  neglected,  using  the  procedures 
described  in  Reference  10.  The  body-fin  junction  leading  edge  treatment 
described  above  is  then  applied.  In  some  cases  the  turning  angle  predicted  by 
(31)  becomes  excessive,  producing  an  unrealistically  large  pressure  jump  and 
possibly  subsonic  flow.  Since  viscous  effects  are  likely  to  be  important 
in  these  regions,  the  computational  algorithm  allows  a  reduction  of  the 
turning  angle  calculated  from  (31).  For  this  purpose  a  multiplying  factor, 
is  introduced,  the  value  of  which  is  user-selected  and  is  <  1 . 
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If  Mj^  is  significantly  less  than  unity,  the  jump  conditions  assigned  by 
option  1  often  do  not  give  satisfactory  results.  This  is  particularly  true  on 
the  expansion  surfaces  where  the  local  analysis  produces  extremely  low 
pressures.  The  condition  M^<< 1  occurs  on  highly  swept  wings  which  feature 
leeside  separation  and  a  large  leeside  vortex.  This  vortex  generates  a  suction 
on  the  upper  wing  surface  which  increases  the  wing  lift.  Under  these 
circumstances,  the  streamlines  on  the  upper  and  the  lower  wing  surfaces  are 
directed  outward  at  locations  near  the  edges.  The  conditions  at  the  leading 
edge  are  thus  strongly  influenced  by  the  flow  on  the  fin  surface.  To  handle 
this  situation,  option  2  is  introduced.  When  an  interior  point 
moves  onto  the  fin,  the  leading  edge  pressure  and  density  is  determined  by 
averaging  values  at  (Xn-i.Ym)  and  (Xn,Y5,).  Using  (7)  the  magnitude  of 
the  velocity  vector  at  Che  leading  edge  is  determined  and  this  vector  is  assumed 
to  be  parallel  to  the  fin  edge.  The  resulting  expressions  for  the  properties  at 
a  leading  edge  point  (Xn.Y,,)  are : 

Pn,m  “  ^Pn-l,m  Pn,m^'^2 

*’n,m  “  (Pn-l,m  ^n,m)'^“ 

q  =yj2(H^  -  h)  (32) 

2  1^2 

qLz/(l  +  Lz) 

0 

2  1/2 

q/(l  +  Lz) 

2.4.2  TRAILING  EDGES.  At  a  trailing  edge  the  two  points  on  Y=Yf 
representing  the  upper  and  lower  fin  surfaces,  are  coalesced  into  a  single 
interior  flow  field  point.  A  local  analysis  is  used  to  determine  the  flow 
downstream  of  the  trailing  edge.  The  computational  algorithm  proceeds  by 
completing  the  step  in  which  a  trailing  edge  is  encountered  without  taking  the 
fin  edge  into  account.  The  resulting  flow  properties  on  the  upper  and  lower  fin 
surfaces  are  the  upstream  values  and  represent  the  flow  properties  on  the  two 
fin  surfaces  Immediately  upstream  of  the  trailing  edge.  The  local  analysis  uses 
these  two  sets  of  flow  properties  in  conjunction  with  the  local  fin  geometry  to 
determine  the  value  of  the  flow  immediately  downstream  of  the  trailing  edge. 

The  trailing  edge  local  analysis  is  dependent  on  the  Mach  number  normal  to 
the  trailing  edge.  If  the  flow  component  normal  to  the  trailing  edge  on  botli 
surfaces  is  sufficiently  supersonic,  the  streamlines  from  the  upper  and  lower 
sides  of  the  fin  will  turn  onto  a  slip  surface  with  normal  s  (see  Figure  3)  by 
means  of  a  system  of  oblique  shocks  and/or  expansions  which  are  attached  to  the 
trailing  edge.  The  orientation  of  the  slip  plane  onto  which  the  streamlines  are 
turned  is  such  that  the  final  pressures  on  both  sides  of  the  slip  surface  will 
be  the  same.  Reference  13  describes  an  iterative  procedure  for  determining  the 
plane  orientation.  Unfortunately  this  scheme  is  cumbersome  to  apply  and 
ccn-ergence  cannot  be  guaranteed  .  Thus,  this  procedure  lias  been  discaided  in 
favor  of  a  simpler  method  that  turns  both  of  the  surface  streamlines  onto  the 
fin  plane  (l.e.,  s  =  ) .  The  coalesced  property  values  for  p,p,u,v 


^'n, m  “ 
'*’n,m 
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are  then  determined  by  averaging  the  results  on  the  upper  and  lower  surface 
streamlines.  The  velocity  component  w  is  solved  for  from  (7)  which  ensures  that 
the  coalesced  trailing  edge  point  has  the  correct  total  enthaply. 

If  the  flow  components  normal  to  the  trailing  edge  on  either  wing  surface 
are  subsonic,  a  different  algorithm  is  applied.  At  a  trailing  edge  point 
(Xu,Yjj),  flow  properties  are  assigned  to  be  those  at  (Xn+l  Yn,) ,  unless  this 
point  is  a  fin  point.  In  that  case  properties  at  (X^-liYm^  are  used.  If 
both  of  these  points  are  on  the  fin,  the  upper  and  lower  surface  properties  are 
set  to  the  average  values  at  points  (Xp,Yj,+j)  and  (Xj,,Yuj_2). 

At  trailing  edge  points  located  on  the  body-fin  junction  the  flow  is  turned 
within  the  body  tangency  plane  using  shocks  or  expansions  until  it  is  directed 
along  the  fin  plane  with  normal  s  =  ei^  .  The  required  turning  angle  is: 


(33) 


On  an  upper  fin  surface  a  compression  is  required  if  v  <  0  and  an  expansion  if 
V  >  0.  The  reverse  sign  convention  applies  on  lower  surfaces.  The  resulting 
values  of  p  and  p  are  determined  by  averaging  the  downstream  properties  on  the 
upper  and  lower  surfaces.  Since  q*^  “0,  v=0  and  u  and  w  components 
are  solved  for  by  satisfying  the  total  enthalpy  condition  (7)  and  the  body 
boundary  condition  (13).  This  yields: 

2(H  -  h) 
w  =  ■  o  _ 

(1  +  bj.^) 


u  =  bzw 

2.4.3  SURFACE  SLOPE  DISCONTINUITIES.  Fin  surfaces  of  design  interest 
often  contain  slope  discontinuities.  Particularly  when  the  surface  velocity 
component  normal  to  the  discontinuity  is  supersonic,  it  is  advantageous  to  apply 
a  local  analysis.  This  procedure  turns  the  flow  across  the  discontinuity  by 
means  of  an  oblique  shock  wave  or  expansion.  The  computational  algorithm  for 
handling  fin  surface  discontinuities  is  analogous  to  that  developed  for  body 
discontinuities  in  Reference  10.  The  scheme  proceeds  by  completing  the  step  in 
which  the  slope  discontinuity  is  detected  without  taking  into  account  the 
presence  of  the  discontinuity.  The  resulting  values  at  the  end  of  this  step  are 
the  upstream  properties  which  in  conjunction  with  the  local  fin  geometry  are 
used  via  the  local  analysis  to  determine  downstream  conditions. 

A  fin  surface  slope  discontinuity  is  defined  to  occur  if  the  following 
inequality  is  satisfied: 
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where 


-  iZ  Mar.(}v||  ,  )  |  +  |  |e"  -  -  iZ  Ma;:(|3||  ,  |e^-^)  j  ^  2(10“ 

(35) 


=  \;^  -  V.'r  5z  =  ^  e;- 

Xr  X 


The  upstream  flow  properties  and  the  vectors  normal  to  the  fin  surface  upstream 
( nf  )  and  downstream  (nj_|_)  of  the  surface  slope  discontinuity  form  the 
nec'essary  inputs  for^the  turning  procedure  outlined  in  the  following  section. 
The  vectors  nf  and  n^  are  illustrated  in  Figure  3. 


2.4.4  EXPANSION  AND  COMPRESSION  TURNS.  As  indicated  in  Sections  2.4.1  to 
2.4.3,  the  treatment  of  leading  edges,  trailing  edges  and  surface  slope 
discontinuities  involves  turning  the  flow  by  either  a  compression  or  expansion. 
This  section  provides  the  analysis  used  in  implementing  these  turns  and 
constructing  the  velocity  downstream  of  the  turn. 


The  parameters  needed  for  calculating  a  compression  or  expansion  turn  are: 
n_,n+,p_,p_  and  q_.  Here  the  minus  and  plus  subscripts  refer  to  the 
quantities  upstream  and  downstream  of  the  turns  respectively.  The  vectors  n_ 
and  nj.,  which  are  Involved  in  this  analysis,  are  illustrated  in  Figure  5.  In 
the  case  of  a  leading  edge,  the  n_  vector  used  in  this  section  corresponds  to 
the  n_  vector  defined  in  Section  2.4.1  while  the  h+  vector  corresponds  to 
nf .  Similarly,  at  the  trailing  edge  n_  corresponds  to  "Sf  and  "fij.  to  ■§ 
while  at  surface  slope  discontinuities  n_  corresponds  to  nf  and  to 
nf+.  Using  this  information  the  following  turn  parameters  'are  calculated: 


|n_*  n+ 
_ 

h_l 

t"  =  X  if 
s  =  n  X  r 


q 


n 


Rn 

~a 


(36) 


Here  i  is  the  turning  angle,  q^  is  the  upstream  velocity  component  tangent 
to  the  edge  of  the  discontinuity  and  qj,  is  the  upstream  velocity  component 
normal  to  the  edge.  The  vectors  t  and  ¥_  are  illustrated  in 
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Figure  5.  To  facilitate  consideration  of  both  upper  and  lower  fin  surfaces,  the 
parameter  fl  is  introduced  which  has  a  value  of  +1  and  -1  on  the  upper  and  lower 
surfaces  respectively.  A  turn  is  an  expansion  if  <  0  and  a  compression 

if  fl  qn  ^  0.  ~ 

For  an  expansion  turn  the  Prandtl-Meyer  relations  are  applied  if  Mn  ^  /l.05. 
These  are  given,  in  differential  form,  by: 

,  -p  q^ 

"  n  -  (37) 

da 


dp  >  _1  ^ 
da  ^2  da 


-  2(Ho  -  h)  -  q2 

Integration  of  (37)  from  a  *  0  where  qn  “  Qn  »  P  “  P_>  P  “  P_  to 
a  =«  4 ,  using  the  second  order  improved  Euler~ method,  gives  p+, 
p+,qn+.  If  0<  Mn  </1.05  the  flow  is  isentropically  expanded  to 
MnVi.05  using  the  relations: 


where 


Sl 


1  +  (Lli)M2 
2  1 

,  2 
1  +  (Lli)  M 
2  2 


Y 

(Y-1) 


P  * 

c 


q*  = 


1/2 


Ml  -  Mj,  and  M2  »  /1.05. 


(38) 


Equations  (37)  are  then  integrated  using  p*,  p*  and  q*  as  the  initial 
conditions.  The  choice  of  Mj,  **  / 1 .05  as  the  lower  bound  for  the  direct 
application  of  the  Prandtl-Meyer  relations  is  not  crucial  and  any  number 
slightly  greater  than  unity,  which  avoids  special  treatment  of  the  improper 
integral  when  Mj,  =  1,  would  be  satisfactory. 


29 


NSWC  TR  81-457 


A  compression  turn  occurs  if  ^  0.  In  these  cases  several  possible 
turning  strategies  are  considered.  The  domain  of  application  of  each  strategy, 
as  a  function  of  normal  Mach  number,  Mj,,  and  turning  angle,  6,  is 
illustrated  in  Figure  6.  The  oblique  shock  relations  are  used  if  there  exists 
an  oblique  shock  solution  with  supersonic  downstream  flow.  Such  a  solution  is 
permitted  if  Mj,  >  1  and  5  <  5*,  where: 


6* 


tan~^ 


M2sln2a*  -  2coto  * 
n 


2  + 


Mn(Y 


+  cos  2o  * ) 


rj9) 


a  *  *  sln~^ 

(Y  +  1)M2  -  3  +  y  +-yjc 

4y 

C  =  (Y  +  1)  j^(Y  +  1)M^  -  2(3  -  y)  m2  +  Y  +  9j  . 

If  an  oblique  shock  solution  exists,  the  shock  angle  is  determined  by  solving 
the  cubic  relation: 

X^  +  cx2  +  cX+c  =  0 
1  2  3 

c  «  -  (1  +  2/m2)  -  Y  sin2e 

1  n 

c  =  (2m2  +  1)/M^  +  [(y  +  1)2/4  +  (y  -  1)/M2]sin2e 

2  n  n  n 


c  =  cos2$/M^ 

3  n 

for  the  middle  root.  Here  x  is  the  sine-squared  of  the  shock  angle.  With  the 
shock  angle  determined,  the  Rankine-Hugonlot  relations  can  be  applied  to  the 
calculated  conditions  downstream  of  the  shock: 

^  =  [2ym2x  -  (y-1)]/(y+1) 

P_  n 


1+  =  (y+1)m2x/[(y-1)M^+2]  (^1) 

p  n  n 


q^  ll  -  4(M  X  -  1)(ym2x  +  1) 
*  =.  I  n _ n _ 

\  2  4 

(Y+l)MnX 
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If  Mn  >  1  and  an  oblique  shock  is  not  permitted,  the  downstream  flow  field 
conditions  are  determined  using  an  empirical  procedure  which  was  developed  using 
the  data  of  Reference  14.  An  effective  shock  angle  Is  estimated  using  the 
following  relation: 

0  -L-xfl-a*]  (42) 

eff  2  12  J 

X  -  Max  (0,Mln[1.94  -  .65  (L-) ,  3.4  -  1.405(1-)]] 

(S  (S  ^ 


When  6/6*  »  2.42,  the  above  equation  predicts  a  normal  shock.  For  further 
increases  In  6/5*,  the  subsonic  flow  behind  the  shock  Is  Isentroplcally 
compressed  from  Mg,  the  Mach  number  behind  the  shock, to  a  final  Mach  number  of  M^ . 
The  quantity  M^  is  determined  from  the  following  equation: 

Mf  -  Max  [Mg[l  -  111  +  (.13)sln(in(.)],0.] 


=  (6  -2.426*) _ 

Max[. 314, Min(. 855,1. 438-.  815M^^)] 


(43) 


Downstream  conditions  are  then  predicted  by  applying  (38)  with  M]^  “  Mg  and 
M2  *  Mf.  When  M^  <  1,  a  shock  Is  not  associated  with  a  compression  turn 
and  downstream  conditions  are  predicted  by  isentroplcally  compressing  the  normal 
flow  component  to  Mf. 


The  condition  Mjj  <  0  can  occur  on  either  an  expansion  or  compression 
turn  and  Indicates  that  the  velocity  vector  does  not  cross  the  edge  of  the 
discontinuity.  This  situation  is  treated  by  isentroplcally  stagnating  the  flow 
component  normal  to  the  leading  edge  using  (38)  with  Mi  =  Mjj  and  M2  “  0. 

At  the  completion  of  this  procedure,  qn+  =  0  and  the  velocity  vector  is 
directed  along  the  edge  of  the  discontinuity. 


At  the  completion  of  the  expansion  or  compression  turning  procedures, 
p+,p+, qjj+  have  been  determined  and  q-j  is  unchanged.  The  downstream 
velocity  vector,  can  be  constructed  from: 


^  .  q  f 

+  T  1"^  i 

|s+l 

where  s+  •  n_^  x  t 

Here  the  upper  sign  is  used  for  a  compression  turn  and  the  lower  for  an 
expansion  turn. 

2.5  TREATMENT  OF  INLETS 


(44) 


Provisions  are  Included  in  the  computer  program  for  considering  Inlets 
which  are  attached  to  the  missile  body.  The  inlet  may  have  an  arbitrary  shape 
provided  that  the  inlet  cowling  and  the  missile  body  together  can  be  represented 
by  r  “  b(i>,z),  where  b(i>,2)  is  single-valued  in  ♦  for  each  z.  Also  the 
cowling  lip  must  be  sharp  and  must  lie  in  a  plane  perpendicular  to  the  body  axis. 
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The  presence  of  an  inlet  Is  checked  at  each  ♦  plane,  at  every  step.  An 
inlet  is  defined  to  occur  on  plane  m  at  step  k  if: 


bk  _  bk-1 


-  max(  I  I,  |b^J-|)4z 


.01 


(A5) 


where 


Y 

bk  =  b'^  -  z  hk 


Here  b  is  the  body  radius,  b^,^  are  the  derivatives  of  b  with  respect  to 

z  and  if.  The  computational  step  in  which  an  inlet  is  detected  is  completed 

using  the  body  radius  from  the  previous  step,  b^~^.  The  flow  field  at  a  radius 
of  less  than  b^  is  then  excluded  from  the  computation  by  rezonlng  which  places 
the  body  surface  of  plane  m  at  a  radius  of  b^  and  leaves  the  shock  location 
unaltered.  The  same  number  of  radial  points  are  retained  along  plane  m  and 

their  distribution  is  determined  by  f(X,Y,Z)  of  (3).  The  flow  properties  at 

the  new  grid  point  locations  are  determined  from  old  values  by  linear 

interpolation.  At  the  wall  the  Interpolated  properties  are  interpreted  as  the 
conditions  immediately  upstream  of  the  Inlet  lip.  The  conditions  immediately 
downstream  of  the  lip  are  obtained  using  the  local  analysis  of  Section  2.4.4. 
The  velocity  vector  is  turned  through  the  angle: 


6 


(46) 


-*■  A  -k  1  “a  ^ 

Here  n_  =  e  -  _z.  e  +  "  -u)®  and  Is  perpendicular  to  the  plane  containing 

r  b  if  w  b  2 

the  lip  tangency  vector  and  the  velocity  vector  q.  To  accomplish  this  turn  the 
upstream  velocity  vector  is  broken  into  components  normal,  qj,  ,  and  parallel 
qx ,  to  the  lip  using  and  rtb,  the  inlet  surface  normals.  TTTe  resultant 
relation  are: 


*  X 

q* (n_x  T ) 

- =v, 

n_  lb_x  T  I 


(47) 


The  properties  P_fP_,qn_  known  and  serve  as  the  initial  conditions  to  the 
turn  which  is  accomplished  using_^the  procedure  described  in  Section  2.4.4.  A 
compression  turn  is  required  if  q*  nb  0  and  an  expansion  if  q*  nb  >  0. 

The  turning  procedure  produces  p+,p+,qxi+  while  leaving  qx  unchanged. 

The  final  velocity  vector  can  Chen  be  constructed  from 
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q 


n  XT 

n_ 


(48) 


in  which  Che  plus  sign  is  used  for  a  compression  turn  and  minus  sign  for 
expansion. 
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SECTION  3 

SPECIAL  PROCEDURES 


The  presence  of  fin  surfaces  in  the  interior  of  the  computational  domain 
requires  the  introduction  of  some  special  differencing  procedures.  In  addition, 
physical  considerations  have  motivated  other  adjustments  to  the  differencing 
used  at  both  fin  and  interior  points  located  next  to  the  fin  edge. 

3.1  ALTERATION  OF  X  DIFFERENCING  FOR  FIN  AND  INTERIOR  POINTS  ADJACENT  TO 
THE  FIN  TIP.  The  types  of  points  under  consideration  in  this  subsection  are  A 
B,C,  of  Figure  2.  Selection  of  appropriate  schemes  for  advancing  these  points 
depends  on  the  Mach  number  normal  to  the  leading  edge  and  the  applied  leading 
edge  treatment.  Several  differencing  strategies  are  available. 

Option  0  for  the  fin  points  such  as  A,B  and  interior  points  such  as  C  is 
similar  to  Che  computational  algorithm  applied  elsewhere  in  Che  flow  field.  Fin 
points  A,B  are  advanced  using  the  usual  fin  surface  point  algorithm.  The 
MacCortnack  scheme  for  advancing  point  C  must  be  modified  since  there  are  two 
adjacent  sets  of  flow  values  (i.e.,  points  A  and  B)  corresponding  to  the  upper 
and  lower  fin  surfaces.  Point  C  is  advanced  in  two  separate  calculations  using 
first  the  lower  fin  surface  values  at  A  and  then  the  upper  fin  surface  values  at 
B.  The  resulting  two  conservation  vectors  are  Chen  averaged. 

A  second  strategy,  option  1,  advances  the  fin  edge  points  A,B  without  using 
the  information  at  point  C,  and  interior  points  such  as  C  without  recourse  to 
the  information  at  points  A,B.  To  advance  point  C  using  this  option,  X 
differences  are  taken  in  the  direction  away  from  the  fin  in  both  the  predictor 
and  corrector  steps.  Using  one-sided  X  differences  to  advance  points  A,B  has 
been  found  to  produce  unsatisfactory  results.  Instead,  X  derivatives  calculated 
from  flow  properties  at  A  and  C  or  B  and  C  are  set  to  zero.  The  computational 
algorithm  has  been  constructed  with  sufficient  generality  to  allow  C  to  be 
advanced  using  A,  B  ,  or  A  and  B. 

The  choice  of  which  option  to  use  in  calculating  X  differences  for  fin  edge 
points  and  adjacent  interior  points,  is  primarily  dictated  by  the  treatment 
employed  at  the  leading  edge.  If  the  leading  edge  is  treated  without  recourse 
to  a  local  analysis  (i.e.,  leading  edge  option  0),  points  A,B,C  are  advanced 
using  X  differencing  option  0  which  allows  points  off  the  fin  to  be  influenced 
by  points  on  the  fin.  When  a  local  analysis  is  used  (i.e.,  leading  edge 
option  1),  the  X  differencing  option  1  is  recommended.  The  rationale  for 
selecting  this  combination  of  options  can  be  seen  by  considering  the  leading 
edge  with  an  attached  shock  wave.  Here  the  Mach  number  normal  to  the  leading 
edge  is  supersonic  and  the  interior  points  on  the  fin  plane,  such  as  C,  should 
not  be  influenced  by  the  presence  of  the  fin.  The  use  of  option  1  practically 
eliminates  all  of  the  upstream  influence  of  the  fin.  The  X  derivatives, 
calculated  on  the  fin  surface,  should  reflect  local  surface  variations.  If 
information  at  point  C  is  used  to  advance  A,B,  the  X  derivatives  will  reflect 
property  variations  across  the  shock,  and  the  calculated  values  will  be  greatly 
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in  excess  of  the  local  values.  The  use  of  option  1  essentially  assumes  that  the 
flow  properties  at  the  tip  are  equal  to  those  at  the  adjacent  fin  point. 

In  the  case  of  a  subsonic  leading  edge,  a  more  complicated  choice  is 
necessary.  Here  compression  surfaces  are  advanced  using  both  X  differencing  and 
leading  edge  options  of  0,  while  the  expansion  surfaces  are  treated  with  X 
differencing  and  leading  edge  options  of  1  and  2.  The  interior  point  adjacent  to 
the  fin  tip,  C,  is  advanced  using  X  differencing  option  0. 

3.2  SUPPRESSIOK  OF  Y  DERIVATIVES  NEAR  THE  LEADING  EDGE.  Figure  7 
illustrates  the  calculated  surface  pressures  on  a  fin  in  uniform  flow.  The 
calculated  results  should  exhibit  a  constant  pressure  downstream  of  the  leading 
edge,  but  in  fact  overshoot  the  leading  edge  value.  The  excessive  pressure 
values  aft  of  the  leading  edge  are  a  consequence  of  the  numerical  procedure  and 
the  error  becomes  more  severe  as  the  magnitude  of  the  pressure  jump  at  the 
leading  edge  increases.  Such  inaccuracies  at  the  leading  edge  can  have  a  strong 
influence  on  the  total  vehicle  aerodynamics.  The  calculated  pressure  overshoot 
at  the  leading  edge  may  be  suppressed  by  damping  the  Y  derivatives  which  occur 
in  (25)  and  (27)  that  advance  the  fin  and  corner  pressures.  Such  a  procedure, 
described  below,  is  automatically  implemented  on  leading  edges  which  feature  a 
pressure  rise. 

Following  the  occurance  of  a  leading  edge  the  Y  derivatives  are  set  to  0 
for  one  step  and  damped  for  an  additional  Kg  steps.  The  damping  is 
accomplished  by  multiplying  the  Y  derivatives  by  a  factor  0  £  D,,  £  1 .  A 
new  Kg  is  calculated  for  each  leading  edge  point  encountered.  The  value  of  Kg 
for  a  leading  edge  point  occuring  at  (Xj,,Y^,Z|^)  is  based  on  the  magnitude 

of 


4£ 


k+1  k+1 

Pn,m-*-l  ”  ^n,m 


^ ♦n ,m+l  ^n,m) 


Here  the  top  and  bottom  signs  apply  to  the  upper  and  lower  fin  surfaces 
respectively,  while  the  k+1  superscript  indicates  that  this  quantity  is 
evaluated  one  step  following  the  occurrence  of  the  leading  edge.  The  algorithm 
for  selecting  Kg  is: 


Ks 


.5655,  4.28  +  .08 


^  )  Max 
4<t>/ 


'1.33 


(49) 


where  ([^]]=  largest  integer  not  exceeding  x. 
The  Damping  factor  D^,  is  calculated  from; 


D 


c 


0 


k  <  (Kg  -  4) 


'l 


(X  -  k) 

s  ■_  ■  (K  -  4)  <  k  <  K 

4  s  s 


(50) 
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Here  k  is  the  number  of  stops  following  the  occurence  of  the  leading  edge. 

At  the  body  fin  junction,  the  algorithm  for  predicting  Kg  at  point 
(Xx.Yn,)  is  given  by 


.35 

A4>  . 


If  k  <  Kg  the  corner  relations  for  advancing  pressure  are  not  used. 

Instead,  the  velocity  vector  is  isentropically  turned,  first  within  the  body 
tangency  plane  and  then  within  the  fin  tangency  plane,  until  it  is  parallel  to 
the  corner. 


3.3  APPLICATION  OF  SMOOTHING  TO  INTERIOR,  BODY,  AND  FIN  POINTS.  In 
computations  featuring  body  separation,  and  on  highly  swept  wings  with  subsonic 
normal  Mach  numbers  at  the  leading  edge,  large  vortex  structures  develop  in  the 
flow  field.  In  such  circumstances  it  is  often  necessary  to  smooth  the 
calculated  flow  field.  This  is  accomplished  by  applying  a  switched  Schuman 
filter^^  with  a  density  switch  after  the  completion  of  each  corrector  step, 
prior  to  decoding.  (The  use  of  the  Schuman  filter  is  one  way  of  introducing 
artificial  viscosity.)  For  interior  points  smoothing  is  applied  to  the  con¬ 
servation  vector  U,  while  at  the  body  and  fin  surfaces  it  is  applied  along  the 
surface  to  the  advanced  quantities.  The  smoothing  algorithm  is  as  follows: 


interior 


(52a) 


Ui,j  =  ti,j  +  ^^i+l,j  ~  ^i,j)Ci+i/2,j  -  (^i,j  -  i^i-l,j) 

(^i,j+l  “  ^,j>Ci,j+i/2  -  (^,j  -  i55,j-l)Ci  j_i/2 


.c 


Wi,j  =  Wi  j  +  (Wi^j+l  -  Wi^j)  Ci^j+i/2  -  (^J^i.j  -  ^i,  j-l^Ci  j_i/2 


*^i,j  °  Hi  j  +  (Hi+ij  -  Hij)Ci+i/2  j  -  (Hij  -  Hi_i  j)Ci_i/2,j 


(52b) 


(52c) 


(52d) 


•^i,j  *  ^i,j  1  ^^i,j±l  “  >^1,  j^^i  j+i/2  -  (J  -  J  i+l,j^^i+l/2,j 


37 


NSKC  TR  81-457 


1 


[  Pi.j  ) 

In  Pi,j  j 

1  •  *  ( 

1 

j  ( 

i  .  1 

u  _li  +  V  1  1 

u  +  (rj  )v  .  .  \ 

'  i.j  bj  i,i' 

'  •  '  V  /  1,1  1 

C 


(pi+i.j  *  pi.j) 


c 

X 


In  Pi,j 


®i.j  ) 

C 

i.j+1/2 


C 

Y 


Here  C^iCy  are  specified  constants.  The  superscript  c  indicates  the 
quantity  values  following  the  corrector  step.  The  switched  Schuman  filter 
smooths  quantities  in  conservative  form  which  preserves  the  numerical  scheme's 
ability  to  capture  shocks  at  interior  points.  Application  of  the  switched 
Schuman  filter  to  the  advanced  wall,  body,  and  fin  quantities,  rather  than  to 
p,p,u,v,w,  ensures  that  the  final  flow  variables  at  a  surface  point  satisfy 
the  flow  tangency  boundary  condition  and  the  total  enthalpy  constraint,  (7), 


Smoothing  is  applied  only  when  requested  by  the  user  or  at  interior  points 
if  the  decoded  pressure  value  is  negative.  In  the  later  case, only  the  point  at 
which  the  negative  pressure  occurred  is  smothed  using  "  *25  and  ®  0. 

A  diagnostic  message  is  also  printed  indicating  that  this  procedure  has  been 
applied.  Ln  the  sample  calculations  discussed  in  Section  4  the  user  specified 
smoothing  is  applied  only  where  indicated  and  the  automatic  smoothing  is  rarely 
invoked. 


3.4  BODY  CROSSFLOW  SEPARATION 


At  incidences  greater  than  about  10°,  the  flow  on  the  leeside  of  a 
slender  body  separates  and  forms  a  symmetric  vortex  pattern.  When  such  a 
configuration  is  treated  with  an  inviscid  code,  a  crossflow  shock  develops 
(illustrated  qualitatively  on  Figure  8)  and  the  surface  pressure  ahead  of  the 
shock  becomes  unrealistically  low.  It  has  been  observed  in  Ref.  16  that  a  more 
realistic  lee  side  flow  pattern  may  be  achieved  by  prescribing  body  surface 
boundary  conditions  that  facilitate  crossflow  separation.  This  section 
describes  a  procedure  that  may  be  used  with  the  present  code  to  simulate 
crossflow  separation  in  cases  where  pitch  plane  symmetry  exists.  The  present 
method  prescribes  the  direction  of  the  two  body  surface  streamlines  located 
nearest  to  the  experimentally  observed  separation  line. 

The  separation  line  on  the  body  surface  can  be  specified  bv  prescribing 

♦  s  =  'Cs  (53) 

The  position  vector,  r g ,  of  the  point  on  the  separation  line  in  the  meridian 
plane  i^g  is  given  by 

Tg  =  b(*s>  ®r  2 
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and  Che  unit  vector,  Xg,  tangent  to  Che  separtion  line,  by 


_  dr  dr 

T  =  _ 1  /  1  _ ® 

s  dz  ■  dz 


(54) 


where 


dr  3  r  ar  _ 

_ ^  ♦ '  +  _ 2.  =  (i>'b  +  b)e  +  4'be  +  e 

dz  3^  s  3z  s  ^  z  r  s  ^  z 


d^ 

♦  '  =  _ 2 

s  dz 


The  specific  form  of  (53),  used  in  the  present  program,  is  the  following 
correlation  of  experimental  data  for  circular  bodies  (Reference  17): 

4  =  2.23  ii 

s  (2 


z/b(z)  -  3. ] tana i  “‘23 


(53a) 


for  z  >  (.450/ tan  a  +  3)b(z). 

Noting  that  the  crossflow  velocity  is  the  projection  of  Che  velocity  vector 
onto  the  z  =>  constant  plane,  and  that  the  tangent  to  the  body  in  the  crossflow 
plane  is  (b^/b)  e^.  +  «(>  >  the  ratio  of  the  total  velocity  to  the 
crossflow  velocity  is  given  by 


[(b^/b)  -f  ^  ^ _ 

qV(b^/b)2  +  1  ■^(2H^  -  h)  [(b^/b)2  +  l] 


(55) 


At  the  separation  point,  Che  surface  velocity  is  prescribed  to  be  in  the 
direction  of  Xg  and  is  given  by  (54).  This  is  enforced  by  solving  (55)  for 
V2  using  1(1  =  i|(g,  where  <(ig  is  determined  from: 


Rb^/b)  1  +  1^]  •  ? 

9  r  9  s 
[(b  /b)2  +  111/2 

♦ 


{0 


/b^/b) 


(h/h)  (9' 

9 _ S 

+  1]  [(^'s 


b^  +  b  )  +  9 '  b 
9  z _ s _ 

2  2  2 

b9  +  bg)  +  9's  b 


+ 


(56) 


since  the  present  method  does  not  explicitly  track  the  separation  point,  (56) 
cannot  be  imposed  directly  because  separation  generally  occurs  between 
computational  points.  It  is,  therefore,  assumed  that  the  crossflow  velocity 
ratio,  9,  varies  smoothly  on  Che  body  surface  and  that  the  separation  affects 
only  the  two  body  points  adjacent  to  the  separation  point. 

Consider  the  case  where  the  separation  point  lies  between  (X|,  Y^,)  and 
(Xi,  The  required  9's  at  (X^,  Y^)  and  (X|,  Y^+i)  are 

determined  using  the  following  linear  interpolations: 
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'•'l.n  "  'i'l.m-l  ^'*'s  "  "  ♦  1  ,m-l  ^ s  ~  '•’l.m-l' 

(57) 

'I’l.in+l  =  ’t’l,m+2  (’^'s  "  *1,01+2^  "  *  1  ,ni+2  s  ~  ■<’l,m+2'^ 

Here  and  i(i  n,+2  are  calculated  by  evaluating  (56)  using  the  flow 

properties  at  points  (X^,  Y^_x)  and  (Xn,  Yni+2 )  respectively.  To  ensure 

that  pressure  and  density  behave  smoothly  near  the  separation  point,  s  and  p  at 

points  (X]^  and  (X]^,  Yjj,+  |)  are  determined  by  linear  interpolation; 

^l.nr+j  °  “l.m-1  (Pl,m+2  “  Pl,m-1^  (♦l.rn+j  *  1  ,m-l )/ (<>  1  ,m+2  "  ♦l.m-l) 

(58) 

®l,ad-j  =  ^l,m-l  (®l,m+2  “  ®l,m-l)  (♦l.m+j  “  -)>  1  .m-l )/ U  l  ,ni+2  '  ♦i.ni-l) 
where:  j  =  1,2 

At  the  end  of  each  computational  step  the  advanced  quantities  at  points  (X]^, 

Yjj,)  and  (Xj^,  Y^^^^)  are  decoded  using  (17)  which  enforces  the  tangency 
boundary  conditions  and  the  correct  stagnation  enthalpy. 

Application  of  (55),  (56),  (57)  and  (58)  on  a  circular  body  produces  a  vertex  i 
the  leeward  flow  field.  Crossflow  velocities  which  are  positive  on  the 

windward  side  of  the  body  become  negative  leeward  of  the  separation  point.  When 

the  flow  field  vortex  becomes  sufficiently  large,  an  interval  of  positive 
surface  crossflow  velocities  (i.e.,  a  secondary  separation)  develops  between  the 
separation  point  and  the  leeward  s^ymmetry  plane.  This  interval  is  unstable  and, 
shortly  after  such  a  region  develops,  the  calculation  fails  unless  the  body 
smoothing  option  is  applied. 
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SECTION  4 

COMPARISON  OF  CALCULATION  AND  EXPERIMENT 


This  section  presents  comparisons  of  calculation  and  experiment  using  the 
algorithms  outlined  in  Sections  2  and  3.  A  description  of  the  computer  code, 
listing  and  set  of  user's  instructions  can  be  found  in  Ref.  11.  The  range  of 
cases  covered  in  this  section  includes  body  alone  with  and  without  crossflow 
separation,  body-fin  configurations  with  supersonic,  subsonic  and  transonic 
leading  edges  and  a  body-wing-tail  configuration  with  and  without  tail 
deflection.  The  input  variables  used  to  execute  each  run  are  listed  in 
Ref.  11.  Most  of  the  cases  presented  in  this  section  have  been  reported 
previous ly^"^  ’ ’  .  However,  different  versions  of  the  numerical  method  were 

used  in  these  references  and  in  some  instances  the  results  may  differ  from  those 
shown  in  this  report. 

All  the  cases  featured  in  this  report  have  bodies  with  pointed  nose  tips. 
The  initial  data  plane  is  located  near  the  tip  and  the  flow  field  on  this  plane 
is  generated  using  a  method  described  in  Ref  11.  This  procedure  provides  an 
estimate  for  the  flow  field  about  a  cone  tangent  to  the  body  at  the  starting 
plane  by  integrating  one-dimensional  conical  equations.  The  resulting  flow 
field  is  approximate  except  at  zero  incidence.  Experience  suggests  that  the 
calculations  described  in  this  section  are  relatively  insensitive  to  the  flow 
field  used  at  the  starting  plane. 

4.1  BODY  ALONE .  A  body  alone  calculation  was  carried  out  with  and  without 
crossflow  separation  on  a  tangent  ogive  cylinder  of  nose  fineness  3,  and 
afterbody  length  of  10  calibers.  The  angle  of  incidence  was  15°  and  the  free 
stream  Mach  number  was  3.  To  minimize  the  computational  cost,  the  runs  were 
completed  in  three  steps.  The  first  part  covered  Z/D  <  2  and  used  a  uniform 
11  X  13  mesh*  while  the  second  and  third  used  a  31  x  37  mesh.  The  third  section 
covered  Z/D  >  5  and  featured  mesh  clustering  in  the  r  direction.  In  the 
crossflow  separation  run,  the  separation  started  at  Z/D  =2.  As  shown  in  Fig. 

9,  application  of  the  crossflow  separation  brings  the  flow  field  into  reasonable 
qualitative  agreement  with  the  measurements  of  Oberkampf^®.  Calculated 
surface  pressures  with  and  without  crossflow  separation  are  compared  in 
Fig.  10.  In  both  cases  the  computed  pressure  is  considerably  lower  than 
experimental  in  the  vicinity  of  separation.  The  pressure  rise  on  the  leeward 
side  of  the  body  in  the  crossflow  separation  case  is  triggered  by  separation, 
while  in  the  other  case  it  is  due  to  the  crossflow  shock  which  develops  on  the 
leeside  of  the  body. 

*The  notation  CN  X  M)  indicates  N  planes  in  the  r  direction  and  M  planes  in  tne 
♦  direction. 
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4.2  CON FIGURiXT ION'S  WITH  SUPERSONIC  LCADIN'G  EDGES.  In  Ref.  21  a 
tangent  ogive  body,  equipped  with  tail  fins  or  several  different  pi  an fc'rms ,  is 
tested  in  supersonic  flow.  The  fins  feature  surface  slope  discontinuities  at 
various  locations  along  the  surface.  The  free  stream  Mach  number  is 
sufficiently  large  to  allow  an  attached  shock  solution  at  the  fin  leading  edge 
in  almost  all  cases.  Numerical  results  have  been  compared  to  experimentally 
measured  surface  pressures  taken  at  Mach  3.7  for  configurations  featuring 
clipped  delta  and  cranked  tail  fins. 

Calculations  for  these  cases  were  made  in  two  sections,  A. uniform  13  X  13 
mesh  was  used  to  advance  the  calculation  to  an  axial  station  slightly  forward  of 
the  tail  leading  edge.  The  flow  field  was  then  rezoned  to  a  uniform  40  X  37 
mesh  for  the  computations  over  the  tail  surfaces.  The  leading  edges  were 
treated  using  option  1  (i.e.,  by  applying  the  local  analysis)  and  X  differencing 
option  1  was  used  until  the  fin  tip  was  encountered.  Along  the  fin  tip  the 
adjacent  interior  point  was  advanced  using  option  0  as  was  the  fin  tip  point  on 
the  compression  (i.e.,  lower)  side  of  the  fin.  The  fin  tip  point  on  the 
expansion  surface  of  the  fin  was  advanced  using  option  1. 

Calculated  and  measured  fin  surface  pressures  are  compared  for  the  cranked 
and  clipped  delta  fin  configurations  in  the  and  "X"  roll  orientation  in 

Figs.  11  to  15.  All  figures  show  reasonable  agreement  between  calculation  and 
experiment.  The  scatter  in  Che  experimental  data  is  a  result  of  using 
experimental  measurements  from  repeated  runs.  On  fin  surfaces  which  feature 
strong  leading  edge  shocks,  the  predicted  results  tend  to  be  lower  chan 
experimental.  Also,  the  pressures  are  not  accurately  predicted  along  the 
body-fin  junction.  The  body  boundary  layer  and  the  fin  leading  edge  shock 
interaction,  presumably  have  an  influence  on  the  corner  and  account  for  much  of 
this  discrepancy.  On  fin  surfaces  which  have  a  weak  leading  edge  shock  or 
expansion,  the  predicted  and  measured' fin  tip  pressure  profiles  are  in  good 
agreement.  However,  the  pressures  along  the  corner  are  underpredicted.  Ch/er 
the  entire  span,  calculated  pressures  on  the  trailing  edge  panel  tend  to  be  less 
than  measured,  probably  reflecting  the  existence  of  a  thick  boundary  layer  or 
separation. 

A  flow  field  vector,  pressure  and  density  contour  plots,  in  a  crossflow’ 
plane  upstream  of  the  trailing  edge,  are  shown  in  Fig.  16,  There  are  no 
experimental  data  available  for  comparison;  however,  the  results  appear  to  be 
qualitatively  correct  and  show  a  leading  edge  shoc'x  propagating  into  the  flow 
field  from  all  three  fins. 

The  normal  force  and  pitching  moment  computations  for  an  airplane  type 
configuration  at  Mach  2  are  shown  on  Fig.  17  in  comparison  to  experimental  data 
of  Ref.  22.  This  case  features  attached  shocks  on  the  wing  and  tail  surfaces. 

The  calculation  was  carried  out  using  a  uniform  13  X  13  mesh  forward  of  the  v^ins 
and  a  uniform  19  X  19  mesh  over  the  remainder  of  the  body.  The  leading  edge  and  X 
differencing  options  were  set  to  1,  To  prevent  the  flow  from  becoming  subsonic 
along  the  wing-body  junction,  the  leading  edge  turning  angle  at  this  location 
x<7as  damped  by  setting  Cp  =  .25  for  wing  compression  and  tail  surfaces  and  =  .5 
on  the  wing  expansion  surface.  The  computed  and  measured  normal  force  and 
pitching  moment  coefficients  for  body-wing  and  body-wing-tail  configurations  at  an 
incidence  of  10*^  are  in  good  agreement.  The  computed  crossflow  field  at  an  axial 
station  just  upstream  of  the  wing  trailins  edge  and  tail  section  are  illustrated 
Fig.  13. 
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4.3  CONFIGURATIONS  WITH  TRANSONIC  LEADING  EDGES.  This  section 
focuses  on  cases  where  the  Mach  number  normal  to  the  leading  edge  is  near  unity 
and  is  too  small  to  allow  an  attached  leading  edge  shock.  Under  these 
conditions  the  local  analysis  at  the  leading  edge  is  approximate. 

The  wind  tunnel  tests  of  Ref.  14  offer  an  opportunity  to  compare 
calculation  and  experiment  for  the  swept  wing  configuration  depicted  in 
Fig.  19.  Measurements  were  taken  on  two  wings  with  the  same  planform  but 
different  thicknesses.  The  calculations  were  performed  at  Mach  numbers  of  2.5 
and  4.5,  at  incidences  of  2°  and  6°  using  a  uniform  15  X  13  mesh,  and  with 
leading  edge  and  X  differencing  options  set  to  1.  The  wing  section  of  the 
configuration  was  also  rerun  on  a  finer  25  X  25  mesh  for  the  thin  wing 
configuration  at  an  incidence  of  6°  at  both  Mach  numbers.  A  comparison  of 
calculated  and  measured  wing  and  body  surface  pressures  is  presented  in  Figs.  19 
to  24.  The  calculated  crossflow  field  of  the  thin  wing  configuration  at  a  = 

6°  is  presented  in  Figs.  25  and  26  for  Mach  numbers  of  4.5  and  2.5 
respectively.  In  each  case  the  crossflow  velocity  vectors,  pressure  and  density 
contours  are  illustrated  at  axial  stations  slightly  forward  of  the  trailing  edge 
point  of  the  wing-body  junction. 

Another  example  of  a  configuration  with  transonic  leading  edges  for  which 
experimental  data  were  available,  is  illustrated  in  Figs.  27  to  30.  In  Ref.  23, 
pressure  measurements  are  provided  on  wing  and  body  surfaces  for  both  and 
"X"  roll  orientations.  Calculations  were  preformed  at  Mach  2.7  at  10® 
incidence  for  both  roll  orientations.  A  uniform  11  X  13  mesh  was  used  forward 
of  Che  wing  and  a  uniform  18  X  21  mesh  for  the  remainder.  Leading  edge  and  X 
differencing  options  were  set  to  1.  Calculated  and  measured  fin  surface 
pressures  are  compared  in  Figs.  27  and  28  which  represent  the  and  "X"  roll 
orientations  respectively.  In  both  cases  the  leading  edge  points  near  the  body 
feature  detached  shocks,  while  those  further  outboard  on  the  wings  have  attached 
shocks.  Figs.  29  and  30  compare  calculated  and  measured  body  surface  pressures 
for  the  "+"  and  "X"  roll  orientations  respectively.  The  calculated  crossflow 
fields  forward  and  aft  of  the  trailing  edge  are  illustrated  in  Figs.  31  and  32 
for  both  of  these  orientations. 

The  final  transonic  leading  edge  example  considered  is  a  swept  wing  model 
depicted  in  Fig.  33.  The  experimental  data  are  taken  from  Ref.  24  and  include 
pressure  measurements  on  the  wing  and  body  surfaces.  Calculations  have  been 
made  at  Mach  numbers  of  2.3  and  2.96  at  incidences  of  8.8°  and  8.6° 
respectively.  A  uniform  13  X  13  mesh  was  applied  to  the  nose  region  (z  <  5) 
while  the  remainder  of  the  body  was  calculated  using  an  18  X  19  mesh.  Option  1 
was  used  for  the  leading  edge  and  for  the  X  differencing.  Wing  surface  pressure 
results  are  compared  to  experimental  values  in  Figs.  33  and  34  for  Mach  numbers 
of  2.3  and  2.96  respectively.  Calculated  and  measured  body  surface  pressures 
are  shown  in  Figs.  35  and  36  and  the  computed  crossflow  plane  flow  fields  at  an 
axial  station  upstream  of  the  trailing  edge  are  illustrated  in  Figs.  37  and  38. 

A  comparison  of  calculated  and  measured  fin  surface  pressures  for  all  three 
configurations  indicates  that  on  expansion  surfaces  agreement  is  reasonably  good 
except  at  the  combination  of  low  Mach  number  and  high  incidence  conditions, 
where  experimental  values  are  smaller  than  the  predicted  ones.  On  compression 
surfaces  best  agreement  was  obtained  for  the  swept  configuration  of  Ref.  14  (see 
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Figs.  19-22).  This  is  to  be  expected  since  these  data  were  used  to  generate  the 
semi-empirical  leading  edge  relation.  Reasonable  agreement  was  obtained  for  the 
other  two  cases  (i.e.,  Figs.  27,  28  and  33,  34)  which  lends  credence  to  the 
semi-empirical  leading  edge  analysis.  The  calculated  body  surface  pressures 
(see  Figs.  23,  24,  29,  30,  35,  36)  also  feature  correct  behavior  in  the  presence 
of  Che  wing.  On  the  windward  side  of  the  body,  surface  pressure  increases 
generated  by  the  presence  of  the  wing,  are  correctly  accounted  for,  although  the 
location  of  the  predicted  pressure  peak,  in  some  cases,  is  downstream  of  the 
measured  one.  On  the  leeside  of  the  body  the  calculated  body  surface  pressure 
decreases  appropriately  in  the  vicinity  of  the  wing,  and  there  is  no  evidence  of 
a  discrepancy  between  the  locations  of  the  predicted  and  measured  pressure 
minimums.  In  some  cases  (i.e..  Fig.  23,  24)  there  is  a  substantial  difference 
between  predicted  and  measured  pressure  levels  on  the  leeside  of  the  body.  This 
is,  most  likely,  attributable  to  the  leeside  boundary  layer  influence. 

4.4  CONFIGURATIONS  WITH  SUBSONIC  LEADING  EDGES.  The  cases  examined  in 
this  section  feature  highly  swept  wings  having  leading  edge  normal  Mach  numbers 
that  are  less  than  unity.  Under  these  conditions,  the  flow  rounding  the  edge  of 
the  wing  separates  and  rolls  up  to  form  a  vortex  which,  in  turn,  induces  a  large 
suction  on  the  lee  surface  of  the  wing.  Previous  studies  (e.g..  Ref.  25)  have 
indicated  that  numerical  solutions  to  Euler's  equations  exhibit  this  type  of 
flew  field  for  wings  with  subsonic  leading  edges. 

The  first  case  considered  is  the  swept  wing  shown  in  Fig.  39.  Experimental 
surface  pressures  are  taken  from  Ref.  26,  and  the  comparison  is  shown  for  a  Mach 
number  of  2.86  and  an  incidence  of  10.3®.  The  current  computational  procedure 
requires  the  presence  of  a  missile  body.  Hence,  it  was  necessary  to  add  to  the 
wing-alone  shape  a  small  circular  centerbody  of  a  diameter  equal  to  the  wing 
thickness.  The  nose  tip  for  this  body  consisted  of  a  cone  of  a  half  angle  of 
5°.  The  presence  of  this  center  body  invalidates  a  surface  pressure 
comparison  near  the  center  of  the  wing.  However,  farther  along  the  span,  this 
wing  should  provide  a  meaningful  comparison  of  experiment  and  calculation.  The 
computation  was  carried  out  with  a  uniform  19  X  19  mesh.  Smoothing  was  applied 
to  the  interior  of  the  flow  field  and  also  to  the  wing  and  leeward  body 
surfaces.  On  the  windward  surface,  option  0  was  used  for  the  leading  edge  and 
the  X  differencing  while  on  the  leeward  surface  the  respective  options  were  2 
and  1.  The  interior  point,  adjacent  to  the  wing  tip,  was  advanced  using  X 
differencing  option  0.  The  surface  slope  discontinuity  jump  was  suppressed  for 
z  <  10.6.  A  comparison  of  the  calculated  and  measured  surface  pressures  on 
the  upper  and  lower  wing  surfaces  is  presented  in  Fig.  39.  The  leeside  pressure 
profile  exhibits  the  expected  suction  near  the  outer  tip.  On  the  compression 
surface  the  predicted  pressures  tend  to  be  lower  than  experimental.  The 
crossflow  velocity  vectors  and  pressure  and  density  contours  are  plotted  on  Fig. 
40  at  an  axial  station  just  upstream  of  the  trailing  edge.  The  expected  leeside 
vortex  is  evident. 

In  Fig.  41  normal  force  and  center  of  pressure  prediction  for  two 
body-wing-tai 1  configurations  are  compared  to  experimental  data  of  Ref.  27. 
Calculations  for  these  bodies  have  been  carried  out  at  a  free  stream  Mach  number 
of  2.86  and  incidences  of  6°  and  ll'’,  with  and  without  tail  deflection.  The 
compi'-tati  ens  vjere  completed  in  three  sections:  upstream  of  the  wing,  wing 
section,  and  the  tail  section.  The  respective  mesh  sizes  were  13  X  13  ,  19  X  25  , 
and  19  X  25.  The  effect  of  the  horizontal  tail  deflection  of  -20°  is  shovTi  in 
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Fig.  42.  The  deflected  tail  case  was  simulated  by  decreasing  the  slope  angle,  v  , 
on  the  compression  surface  by  tan(20“)  and  increasing  it  on  the  expansion  surface 
by  the  same  amount.  The  wing  featured  subsonic  leading  edges,  and  leading  edge 
option  0  was  applied  on  compression  surfaces  and  option  2  on  expansion  surfaces. 

The  X  differencing  option  0  was  used  on  the  compression  surface  and  on  the  interior 
point  adjacent  to  the  fin  tip,  while  option  1  was  used  on  the  expansion  surface. 

The  wing-body  junction  turn  was  damped  by  setting  C.^  =  .5  for  both  the  wing  and 
tail  surfaces.  For  the  section  of  the  calculation  which  covered  the  wing  surfaces, 
interior  points  were  smoothed  using  Cjj  =  Cy  =  .2  while  the  expansion  wing  surface 
and  body  surface  leeward  of  the  wing  was  smoothed  using  j+l/2  “  section 

3  of  the  calculation  which  covered  the  tail,  leading  edge  and  X  differencing  options 
were  set  to  1  and  the  interior  flow  field  was  smoothed  using  Cv  =  Cy  “  *3.  The 
calculated  normal  force  and  center  of  pressure  are  in  reasonable  agreement  with 
experiment  for  both  deflected  and  undeflected  tail  surfaces.  The  crossflow  velocity 
vectors  are  illustrated  in  Fig.  43  for  crossflow  planes  aft  of  the  wing  trailing 
edge,  forwards  of  the  tail  leading  edge  and  aft  of  the  tail  trailing  edge. 
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SECTION  5 

CONCLUDING  REMARKS 


A  numerical  method  has  been  developed  which  predicts  the  inviscid 
supersonic  flow  field  about  finned  configurations  of  engineering  interest.  The 
computational  requirements  are  generally  modest.  All  of  the  cases  demonstrated 
in  this  report  required  less  than  two  minutes  on  a  CDC  7600  and  used  between 
110k  and  170k  core  storage  (octal).  The  present  study  differs  from  previous 
methods  by  treating  the  fin  and  body  geometries  separately.  Simple 
transformations,  based  on  the  body  alone  geometry,  are  used  to  map  the  physical 
plane  into  the  computational  one.  A  local  analysis  is  applied  to  fin  edge  which 
facilitates  the  use  of  coarse  computational  grids.  At  present  a  thin  fin 
approximation  is  employed  which  limits  the  applicability  of  the  computational 
procedure  to  relatively  thin  fins  with  sharp  leading  edges.  With  this 
formulation  it  is  possible  to  treat  a  wide  variety  of  configurations  of 
engineering  interest  which  may  contain  an  arbitrary  number  of  fins  at  small 
angles  of  deflection,  camber  or  dihedral.  By  appropriate  modeling  at  wing  tips 
and  at  estimated  body  separation  points,  it  appears  feasible  to  simulate  flow 
field  vortices. 

Special  procedures  are  used  to  analyze  the  flow  on  leading  and  trailing 
edges.  This  treatment  is  dependent  on  the  Mach  number  of  the  flow  normal  to  the 
leading  edge.  The  treatment  of  the  supersonic  case  is  straight  forward  since 
the  local  analysis  is  exact.  In  the  subsonic  and  transonic  cases  the  local 
analysis  is  empirical  and  several  options  are  available  for  treating  fin  leading 
edges.  Here,  experience  with  similar  cases  for  which  there  is  experimental  data 
may  improve  results. 

The  computational  procedure  has  been  applied  to  a  number  of  different 
configurations.  Good  agreement  has  been  obtained  between  experiment  and 
calculation,  particularly  with  respect  to  aerodynamic  coefficients.  The 
calculated  pressure  is  generally  in  reasonable  agreement  with  experiment.  Here 
deficiencies  are  more  likely  to  occur  near  body-fin  junctions  and  close  to  the 
leading  edges  on  wings  with  detached  shocks.  The  pressure  gradients  on  the  wing 
surfaces  reflect  the  correct  variation  with  leading  edge  Mach  number,  even  on 
the  leeward  side  of  the  wing  in  subsonic  flow.  Here  the  appropriate  suction 
appears  on  the  fin  surface  near  to  the  tip.  Vector  and  pressure  contour  plots 
of  the  flow  field  have  also  been  presented  for  a  number  of  cases.  Although 
there  is  no  experimental  information  for  direct  comparison,  Che  flow  field 
exhibits  the  expected  structures.  For  configurations  with  supersonic  leading 
edges,  shock  waves  can  be  seen  propagating  off  the  edges  into  the  flow  field, 
while  in  the  subsonic  case,  a  large  leeside  vortex  is  seen  to  develop.  In 
several  cases  normal  force  and  pitching  moment  coefficients  have  been  calculated 
and  compared  to  experiment.  Good  agreement  has  been  obtained  even  for 
configurations  with  deflected  tails. 
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Figure  4.  The  effect  of  changing  the  leading  edge  pressure  and  streamline 
direction  on  a  swept  wing  surface  pressure  profile 
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Figure  12.  Calculated  and  measured  fin  surface  pressures  on 
delta  fin  configuration  in  the  "+"  roll  position, 
M  =  3.7,  Ref.  21.  (Zero  reference  shifted  by  1.0 
successive  curve.) 


Figure  18.  Continued 

b)  Tills  station  is  •ipstream  of  the  tall. 
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Figure  21 


O  a  A  O  EXPERIMENTAL 
— . — I-  COMPUTED,  COARSE  MESH 
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.  Calculated  and  measured  fin  surface  pressures  on  a  thin 
wing  configuration  of  Ref.  lA  at  M  =  2.5,  a  =  2°  and  6® 
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COMPUTED,  COARSE  MESH 


Calculated  and  measured  fin  surface  pressures  on  a  thick  swept 
wing  configuration  of  Ref.  14  at  M  =  2 . 5 ,  a  =  2 ®  and  6° 


Figure  25.  Calculated  crossflow  field  on  the  thin  swept  wing  configuration 
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FLgiire  26.  Calculated  crossflow  field  on  tlie  thin  swept  wing  configuration 


Figure  27.  Calculated  and  measured  wing  surface  pressures  on  the  delta 
wing  configuration  of  Ref.  23,  in  the  roll  position  at 
M  =  2.7,  a  =  10°.  (Zero  reference  shifted  by  1.0  for  each 
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Figure  il.  Oaiculated  crossflow  field  on  the  delta  wing  configuration 

of  Ref.  23,  in  the  roll  position  at  M  =  2.7,  a  =  10*^,  Z  =  36.5 


rijjiirc'  32.  Calculated  crossflow  field  on  the  delta  wing  configuration  of 
Ref.  23,  In  the  "X"  roll  position  at  M  =  2.7,  a  -  10^^,  Z  =  36. 


DOUBLE  WEDGE  AIRFOIL 


Figure  33.  Calculated  and  measured  wing  surface  pressures  on  the  swept 
wing  configuration  of  Ref.  24  at  M  =  2.3  and  a  =  8.8*’. 

(Zero  reference  shifted  by  1.0  for  each  successive  curve.) 
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Figure  37.  Calculated  crossflow  plane  flow  field  on  the  delta  wing 
configuration  of  Kef.  24  at  M  =  2.3  and  a  =  8.8”,  Z  =  3( 


Fij;ure  38.  Cn  1  ciliated  crossflow  plane  flow  field  on  tlie  delta  wing 
configuration  of  Ref.  2A  at  M  =  2.96  and  a  =  8.6^^,  /,  - 
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Klpiirf  40.  Calculated  crossflow  plane  flow  field  on  tlie  delta  winp 
configuration  of  Ref.  26  at  M  =  2.86  and  a  =  10.3”.  The 
illustrated  station  is  upstream  of  the  trailing  edge. 
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Figure  43.  Calculated  crossflow  plane  flow  field  for  the  W1  body-wing- tall 
configuration  of  Ref.  27  at  M  =  2.86  and  a  =  12°.  The  forward 
plane  is  downstream  of  the  wing  trailing  edge,  the  middle 
station  is  upstream  of  the  tall  leading  edge,  the  rear  station 
is  downstream  of  the  tall  trailing  edge 
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nomenclature  (Cone.) 


6 

S* 

T 

9 


V 


P 

o(r  ,2) 
a* 


Subscripts 


turning  angle  of  normal  flow  component 
maximum  possible  turning  angle  with  an  attached  shock 
vector  tangent  to  fin  leading  edge  (see  Figure  5) 
angle  between  the  fin  surface  tangent  plane  and  the  fin 
plane  in  the  r  direction 

angle  between  the  fin  surface  tangent  plane  and  the  fin 

plane  in  the  2  direction 

density 

fin  surface  function  (see  Figure  2) 
maximum  attached  shock  angle 

angular  orientation  of  fin  plane  (see  Figure  2) 
angle  between  adjacent  symmjnetry  planes  (if  such  planes 
exist ) 

for  symmetric  problems  and  2r  otherwise 


®>  free  stream  conditions 

downstream  of  an  edge  or  surface  discontinuity 
(see  Figure  5) 

+  upstream  of  an  edge  or  surface  discontinuity 

(see  Figure  5) 

I 

w  wall 

f  fin 
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NOMENCLATURE  (Cent.) 


rif  fin  surface  normal  vector  upstream  of  a 

discontinuity 

n  vectors  used  in  calculating  compression  or  expansion 

jumps  (see  Figure  5) 


(r  ,z) 


number  of  r  planes 
?-n  p 

pressure 

velocity  vector 

velocity  component  normal  to  a  discontinuity 
velocity  component  tangent  to  a  discontinuity 
cylindrical  coordinates  (see  Figure  1) 
entropy 


s 


slip  plane  normal  vector  (see  Figure  3) 


T 

I u ,v ,w) 
(X,Y,Z) 

V, 

"3 

a 

Y 


edge  direction  (see  Figure  3) 

velocity  components  (see  Figure  1) 

compul:ational  coordinates 

u(bj/b)  +  V 

v(ra  )  +  u 
r 

angle  of  attack 
ratio  of  specific  heats 


0 


0' 
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APPENDIX  A 

CHARACTERISTIC  ANALYSIS  OF  THE  FIN  SURFACE 

This  Appendix  presents  a  derivation  of  the  equations  used  for  advancing  the 
fin  surface  points  and  the  point  on  the  fin-body  junction.  The  analysis  does 
not  utilize  the  thin  fin  approximation  and  is  valid  for  any  gas  in  thermodynamic 
equilibrium. 

The  upper  and  lower  fin  surfaces  are  each  independently  described  by 
41  =  o(r,z)  +  $£.  On  the  upper  surface,  Che  flow  region  lies  in  ♦  ^o(r,z)  + 
while,  on  the  lower  surface,  the  flow  region  lies  in  ♦  ^a(r,z)  On 

each  fin  surface  the  inviscid  boundary  condition  (22)  must  be  satisfied.  This 
constraint,  in  conjunction  with  the  full  system  of  equations  (9)  overspecifies 
Che  problem.  The  proper  number  of  independent  equations  necessary  for  advancing 
the  solution  along  Che  fin  surfaces  is  obtained  by  combining  the  admissible 
characteristic  compatibility  conditions  associated  with  (9),  with  (22)  on  the 
fin  surface. 

The  two  families  of  characteristics  associated  with  (9)  are  stream  surfaces 
and  Mach  surfaces.  Each  of  these  families  provide  an  infinite  number  of 
surfaces  on  which  the  characteristic  relations  can  be  written.  Except  when 
considering  points  on  the  body-fin  junction  attention  is  restricted  to  planes 
shown  in  Figure  A-1 .  The  selected  stream  surface  plane  coincides  with  the  fin 
surface  and  the  chosen  Mach  surfaces  Intersect  the  fin  surface  along  a  constant 
2  line.  Two  independent  characteristic  compatibility  conditions  are  satisfied 
on  the  stream  surface.  These  relations  represent  the  propagation  of 
disturbances  for  Increasing  z  along  the  fin  surface.  Since  the  flow  domain  lies 
on  just  one  side  of  each  fin  surface,  only  the  compatibility  condition  for  the 
Mach  surface  lying  within  the  flow  domain  for  z  <  Zg  is  consioered 
admissible.  This  surface  represents  the  propagation  of  disturbances  from  the 
interior  points  to  the  fin  boundary  for  increasing  z.  The  other  Mach  surface, 
which  lies  outside  of  the  flow  domain  for  z  <  Zg,  is  disregarded  and 
replaced  by  the  boundary  condition  (22).  Application  of  admissible 
characteristic  compatability  conditions  at  the  boundaries  in  finite  difference 
calculations  was  first  suggested  by  Kentzer^^-.  This  approach  was  used  in 
Reference  10  to  develop  the  equations  for  advancing  body  surface  and  shock 
surface  points  (17)  and  (19). 

For  the  characteristic  anal^^is  of  a  fin  surface,  it  is  convenient  to 
temporarily  Introduce  an  additioq|pl  transformation; 

t  =  Z  =  z 

n  =  X  =  X(  r  ,11  ,  z)  (A.l) 

C  =  5(X,Y,2)  =  t  _  a  -  .Of 


A-1 
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Note  that  in  the  n  ,?  ,C  coordinates,  n  =  0  corresponds  to  the  body 
alone  surface  and  ?  =  0  corresponds  to  the  fin  surface.  A  quasi-linear  system 
equivalent  to  (22)  is  obtained  by  expanding  (5)  in  terms  of  the  dependent 
variables 

Q  =  (P.u.v.w)*^, 

introducing  the  change  of  variables  (A.l),  and  left  multiplying  the  result  by 
the  non-singular  matrix i 


The  result  is  given  by 


L(Q)  =  Ai2  +  (BiQ  +  (E-F)  =  0 

3n  35  r 


(A. 2) 


where  * 


Ik  riii 

L^qJ 

t  = 


,  »  =<£  (n-  +  Hr  F— 1  +  i  h*  F— 1 

I  3QJ  UqJ  r  ^  [BQJ 

(Z)  {  5  F— 1  +  5  F— 1  +  F— 1  ^  > 

z [3 Q J  L ® ^ J  r  ^  I  3 Q J 

are  the  Jacobian  matrices  of  and 


respectively , taken  with  respect  to  the  components  of  Q  and  subject  to  the  energy 
equation  (7).  Explicit  expressions  for  the  matrices  A,  iB,  and  can  be  obtained 
directly  from  (A. 4),  given  below.  Note  that  the  terms  etc.  are  given 
using  (A.l)  by 


Hg  “  X2>  hr  Xj., 

5  z  =  -<='z»  ^r  =  -^r’ 


(A. 3) 


The  pertinent  facts  concerning  the  theory  of  characteristics  associated  with 
system  of  the  type  (A. 2)  will  be  briefly  reviewed  here;  for  a  more  detailed 
explanation  see  Reference  28,  pp.  577-599.  The  characteristic  matrix  associated 
with  (A. 2)  is 


/A*  (Xi,  X2,  X3)  =  X3/A+X2(B  +  X3  ? 


(A. 4) 


*  In  terms  of  the  transformation  r,  (j),  t 


U  =  UrJ--;  F  =  rJ-l(n2  U  +  n^F  +  (2^)0);  G  =  rJ~t  {fjj  +  r^F  +  (E^^/r)G) 

'  r 

E  =  EJ”^ ;  J  =  n^5({;  “  4rG^ 
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«  (l/a2  -  K/p) 


where 


pA2  -  uaK 


0ti2/T  -  vCK 


pA-i  -  wain 


a  "  Aj^w  +  A2U  H-A^v/r,  A^  '*‘*'2^2  ^3» 

A2  =  nrX2  +5r  *3  =  n^X  2  +  5^X3, 

and  pK  =  (3p/3h)p,  (K  =•  -1/h  for  a  perfect  gas). 

A  surface  't»  (5,n,C)  =  0  is  characteristic  at  a  point  if  its  normal 
at  the  point  satisfies  the  characteristic  condition 


^(^1.  ^2.  ^3)  =  det  |^A*(Xi,  X2,  X3)j=  0 


where  Xi  ■  ii.  X2  ■  1]!L  ,  and  X3  =  liL.  The  characteristic  conoid 
3C  3n  3? 

with  vertex  0  =•  (Cq,  Hq*  Cq)  is  the  envelope  of  all  characteristic 
surfaces  through  0.  The  surface  of  the  characteristic  conoid  is  generated  by 
curves,  called  rays  or  bicharacteristlcs  which  are  the  lines  of  contact  between 
the  characteristic  surfaces  and  the  conoid  they  envelope.  These  curves,  or 
rays,  are  given  by  the  ordinary  differential  equations 


^  ^  dJH  ^  =  IH 

dS  3X]^’  dS  3X2’  dS  3X3 


(A. 6) 


where  S  is  a  parameter.  Each  ray  through  0  is  determined  by  selecting  real 
values  for  X2  and  X3  and  determining  X 3  by  satisfying  the  characteristic 
condition  (A.4)  at  0.  The  characteristic  condition  (A. 4)  for  the  system  of 
equations  (A. 2)  is  given  by 

H  -  Hi  (  Xi,  X2,  X3)  H2  (  Xi,  X2,  X3)  =  0 


where 

3„2 

^2.  ^3)  *  ^-4-  and 

„  2  2  2  2  2  2 
"2(^1»  ^2»  ^3)  -(Aj^  +  A2+A  3/r  )  a 

The  ray  cone  therefore  has  two  sheets,  one  corresponding  to  H^  a  0  and  one 
corresponding  to  H2  =  0.  The  rays  generating  the  sheet  corresponding  to 
Hj^  =  0  are  given,  using  (A. 5),  by 

-  A/w,  ii  -  B/w 
dC  d? 


(A. 7) 
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-►3^  -V 

where  3K  /3-v  =  «.  3W  /3A  =  A  H  grad  n*q,  - i  *  B  S  grad  5*q  . 

1  1  12  3X 

3 

Hence  the  sheet  corresponding  to  =  0  is  a  degenerate  cone  consisting  of  a 
single  ray  through  0.  This  ray  corresponds  in  the  physical  space  (z,<ti,r)  to  a 
streamline*  The  sheet  corresponding  to  f(2  =  0  is  a  true  cone  which 
corresponds  in  physical  space  to  the  Mach  cone.  For  each  characteristic 
surface,  the  lines  of  contact  with  the  cone  (bi characteristic  rays)  have  slopes 
given  by: 

dn  =  (aH  /3X  )/(3H  /ax  ),  ^  =  (3H  /ax  )/(3H  /sx  )  . 

dT  2  2"  2  1  ’  dC  23  21  (A. 8) 


where 


aH  /ax  =  2(ow  -  A  ) 
2  1  1 


V^^2  °  ^  [^(8>^ad  n*q)  -  a2  (A^n^  ^2  '’r'*’  ~2  ^3^'^'^] 

^73X3  =  2  ^cr(grad  ?*q)  -  (A^C^  ^2  ^r  ~2 


(A. 8a) 


On  each  characteristic  surface,  the  system  (A. 2)  reduces  to  one  or  more  scalar 
compatibllltv  conditions  given  by 


I  •  L(Q)  =  0 


(A. 9) 


where  i  is  a  left  null  vector  of  A*  associated  with  the  characteristic 
condition  defining  the  characteristic  surface.  For  the  system  (A. 2),  the  left 
null  vectors  can  be  obtained  from  (A. 4)  by  inspection.  Corresponding  to  =  0 
(i.e.,  0  =  0)  there  are  two  independent  left  null  vectors 


=  (0,w,u,v) 


t  =  (0, 0,  _ 1  ,  -A  ) 

2  r  2 


(A. 10) 


Corresponding  to  ^2  “  left  null  vector  is  given  by 

1 2  =  (per,  wKcr  -  Aj^p,  uKcr  -  A2P,  vKO  ~  pA3/r) 

To  determine  the  left  null  vectors,  bicharacteristic  rays,  etc.,  values  of 
X^,  X2,  X3  must  be  specified.  The  characteristic  conditions  =  H2  =  0 
provide,  for  each  family,  one  relation  between  X^,  X2,  X3.  Since 
both  and  ^2  are  homogeneous  of  degree  2  in  Xj^,  X2  and  X^,  it 
follows  that  =  0  and  H  2  =  0  each  describes  a  one  parameter  families  of 
characteristic  surfaces. 


(A. 11) 
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As  previously  indicated,  attention  is  restricted  to  characteristic  surfaces 
which  are  either  tangent  to  the  fin  surface  or  intersect  the  fin  surface  along 
C  =  Hence  ^2  is  set  to  zero  and  since  v^a>o  the  values  of 

(Xj^/X^)  on  5  =0,  which  satisfy  Eq.  (A. 4),  are  given  by 

»  -  ^  (multiplicity  of  2)  (H  “0) 

w  1 

o 


(C^  +  a6)  X 

2  2 
w  -  a 


(A.12) 


where 


^w^/jgrad  ~  a^C? 


and  X^  is  any  non  zero  value.  Since  on  the  fin  surface  B  *  0,  the 
characteristic  surface  corresponding  to  ■  0  with  X  2  *  0  is  tangent  to 
the  fin  surface.  Since  these  characteristic  surfaces  do  not  leave  the  flow 
domain,  the  associated  compatability  conditions  *  ^(Q)  “  0  and  ^2 
*  L(Q)  =  0  (with  X2  =  0)  are  admissible  on  5*0.  The  final  form  of  the 
associated  compatibility  equations  used  in  the  present  work  are  given  by  (23) 
(corresponding  to  (24)  (corresponding  to  ^l^).  These  can  be 

obtained  by  direct  substitution  of  (A. 10)  into  (A. 9)  with  X2  “  0 
considerable  manipulation  using  the  boundary  condition  (20) .  In  (23) ,  the  entropy 
is  introduced  using  the  thermodynamic  relation 

dp  ~  Pdh  *  -  pTds 

where  T  is  the  temperature.  Note  that  the  derivative —L  and  —  ,  appearing 

3  Z  3X 

in  these  equations,  are  the  same  as  _L  and  —  ,  respectively,  on  5=0. 

H  3n 


For  X2  =  0,  there  are  two  distinct  characteristic  surfaces  (and  left 
null  vectors)  assocltated  with  ^2  =  0  corresponding  to  the  values  of 
X1/X3  given  by  (A.12).  To  determine  the  appropriate  choice  of  sign  in 
(A.12),  the  slopes  of  the  associated  bicharacteristic  rays  through  a  point 
(C,n,0)  are  considered.  Using  (A. 8)  and  (A.12)  yields 

dL  ^  d£/ds  ,  ) 

dC  dt/ds  _a2) 

Since  w>a>o  andi3>wj5  2|,  it  follows  that  the  bicharacteristics  are  as 
indicated  in  Figure  A.l.  For  an  upper  fin  surface,  the  flow  region  lies  in 
5>0  and,  therefore,  the  compatibility  condition  corresponding  to  the  upper 
sign  in  (A.12)  is  the  appropriate  choice.  Analogously,  for  a  lower  fin  surface 
the  lower  sign  is  the  appropriate  choice.  For  both  fin  surfaces,  the 
compatibility  condition  (A. 9)  on  5  =  0  with  ^=^3  is  given  by  (24).  This 
is  obtained,  after  considerable  manipulation,  by  evaluating  (A. 9)  on  5  =0, 
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using  the  boundary  condition  (22)  with  X 2  “  0,  X3  /  0,  and  X j/X 3 
defined  by  (A. 12).  Also,  _L  is  replaced  by  3/3n  *  Yr  (3/3Y) 

3n 

where  Y5  =  Y^iXf/CXj.  +ar  X);cf.,  (A.l)  and  (1). 

The  above  analysis  does  not  apply  to  the  points  (C  ,0,0)  along  the  fin-body 
junction.  These  point  are  special  in  that  they  are  boundary  corners  where  the 
flow  domain  lies  in  the  sector  between  the  body  surface  and  the  fin  upper  (or 
lower  surface:  viz.,  the  quarter  plane  (n  ^  0,  5  ^  0}  or  (n  ^  0,  C  ^  0}). 

The  fin-body  junction  is  modeled  in  a  heuristic  manner  which  leads  to  a  separate 
set  of  equations  for  these  points.  Along  the  junction,  both  (22)  and  (13)  are 
satisfied  which  implies  that  the  junction  is  a  streamline.  It  is  also  assumed 
that  the  junction  is  not  a  vortical  singularity  (i.e.,  the  entropy,  pressure, 
etc.  are  single  valued  along  the  junction).  Under  these  circumstances,  the 
compatibility  condition  (A. 9)  with  ^  implies  that  s  is  constant  along 

the  junction.  Therefore,  because  of  (22)  and  (13),  only  one  other  independent 
relation  is  needed  to  determine  all  the  flow  variables  at  the  junction. 


Of  the  two  remaining  families  of  compatibility  conditions,  only  those 
corresponding  to  ^3  (i.e.,  associated  with  H2  =  0)  are  considered.  From 
these  compatabillty  conditions  an  equation  can  be  obtained  for  advancing  the 
pressure  p  along  the  junction  which  is  independent  of  (22),  (13)  and  s  = 
constant.  On  the_^other  hand,  it  can  be  shown  that  the  compatibility  conditions 
corresponding  to  I2  choice  of  X2  and  X3  (with  X'^2  ^3  ^  0)  is 

not  Independent  of  (22)  and  (13)  when  b2  =  =  0  (a  case  that  we  do  not 

want  to  exclude).  The  specific  choice  of  a  compatibility  condition  from  the 
H2  =  0  family,  is  dictated  by  the  fact  that  the  flow  domain  is  a  sector  and 
only  those  compatibility  conditions  associated  with  bicharacteristic  rays 
through  (5,0,0)  which  lie  in  the  appropriate  sector  for  decreasing  5  should 
be  considered.  Among  the  bicharacteristics  associated  with  H2  =  0  satisfying 
this  condition  there  appears  to  be  no  particularly  convenient  choice.  One 
posslblity  is  the  bicharacteristic  lying  on  the  fin  surface  at  (5,0,0);  i.e., 
35/35  =  0  at  (5,0,0).  From  (A. 8),  this  bicharacteristic  satisfies 
3H2/3X3  =  0.  Simultaneous  solution  of  H2  =  0  and  3H2/3X3  =  0, 
using  (22)  and  (13),  yields  two  possibilities  given  by 


* 

Xi+  =  a  X2f  a(Rn2  -  T5z)  + /T/3  ]/« 


(A. 13) 


2  2’  * 
X3+=-X2[(w  -  a)T+a~T)^Kz±a^^Q//R]/n 


where  d  =  (SR  -  T^)  (w2  -  a^)  +  (R  n  2  -  215  n  +  S5  2), 

z  z  z  z 

R  =  I  grad  C  |  ^ ,  S  =  I  grad  n  I  ^ ,  T  =  grad  n*  grad  5  ,  =  R(w^  -  a^)  +  5  ^  a“ , 

z 

and  X2  ^  0.  To  determine  the  appropriate  choice  to  sign  in  (A. 13),  we 
consider  the  shape  of  the  blcharacteristlc  rays  determined  by  (A. 13).  It 
follows  from  (A. 8)  and  (A. 13)  using  (22)  and  (13),  that  at  (5,0,0) 

* 

dTi/d5  =  -a  [a(Rtiz  -  T5  z)  +  /R  /3  ] /«  . 


(A. 14) 


A-6 
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Since  w>a>o  and  R^T^,  It  follows  that  >  a|Rn2  ”  z  I 

and  thus  the  bicharacteristics  corresponding  to  (A. 13)  are  as  indicated  in 

Figure  A-2.  Therefore  for  both  the  upper  and  lower  fin  surfaces  the  upper  sign 

in  (A. 13)  is  the  appropriate  choice.  The  form  of  the  compatibility  condition 

(A. 9)  on  (Cq,0,0)  with  A  »  A3  and  X3  defined  by  (A. 13)  is  independent 

of  X2  and  given  by  (24a).  Another  possibility  is  to  consider  the 

compatibility  condition  associated  with  blcharacterlstlc  ray  lying  on  the  body 

surface,  n  ■  0;  l.e.,  3t /3C  *  0  at  (5,0,0).  These  blcharacterlstlcs 

satisfy  3H2/3X2  H2  ■  0  which  give  X^  and  X2  in  terms  of  an  arbitrary  X3  /  0 

expressed  by  (A. 13)  with  X2  and  X3,  and  n2,  and  R  and  S  Interchanged.  To  determine 

the  sign,  d5/d§  which  is  given  by  the  right  hand  side  of  (A. 14)  is 

considered  with  the  above  Interchanges.  The  corresponding  bicharacteristics  are 

indicated  in  Fig.  A-2.  Hence  the  appropriate  choice  of  sign  is  the  upper  and 

lower  sign  for  the  upper  and  lower  surface  junctions,  respectively.  This 

compatibility  condition,  evaluated  at  the  Junction,  is  given  by  (24b). 


MACH  SURFACE  CHARATERISTICS  INTERCEDING  FIN  SURFACE  ALONG 
=  CONST  ANT(X,-0) 
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APPENDIX  B 

CALCULATION  OF  AERODYNAMIC  COEFFICIENTS 

The  aerodynamic  forces  acting  on  a  missile  configuration  are  determined  by 
numerically  integrating  the  computed  fin  and  body  surface  pressure  distributions. 
The  sign  conventions  for  the  components  of  the  aerodynamic  force  and  moment  are 
illustrated  in  Fig.  B-1.  The  force  and  moments  are  computed  assuming  that  the 
base  pressure  is  equal  to  pao  and  that  the  moments  are  taken  about  the  point 
z  =  Z(;  on  the  z  axis,  as  is  shown  in  Fig.  B-1.  The  derivatives  with  respect 
to  z  of  the  force  components  are  given  by: 


BFz 

3z 


3Fn 

3z 


Ml 

3z 


3Mz 

3z 


Nf 

A 

Ns 

E 

t 

Fb  bb^  ‘I't)  +  J3  Sn  J  Pb( 

i=l 

n'l 

Nf 

Ns 

E 

1=1 

f 

Pb(b  cos<ti  +  bx  sin(t))d(i)  -  23 

n=l 

Nf 

A 

Ns 

E 

i=l 

1 

p^(b^cos(j)  -  b  sin(|))d<J)  +23 

n=l 

Nf 

-E 

i=l 

•'a 

Pb 

(B-1) 


p^[cos(l)^  -  (raj.)sin(j)^]dr 


(B-3) 


-  “  (raj.)sind'^]+sin(*;g[sin4)£  +  (ra^)cos(})f ]  }rdr 

''in  ^  "  (B-4) 


n=l 


3Mx 

3z 


,  .  3Fy  , 


^  J  p^  b^b^  sln(}!d(l)+ 23  J  Pb  sin(J)g(ra2)rdr  (B-5) 

i=l  -fa  n=l  ri„ 


B-1 
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He.E 

3z  1=1  *<3 


cos(J;  d"}) 


^  >On 

-  E=n/ 


Pbcos6g(rcr2)rdr  (b-6) 


where  p,  =  p  -  p„,  tj)  =  (j)  on  fin  surface. 


of  fin  plane 


The  first  terra  in  each  case  represents  the  body  contribution  while  the  second, 
that  of  the  fins.  In  the  above  equations,  Ng  is  the  nuraber  of  fin  surfaces  and 
Nf-1  is  the  nuraber  of  fins  not  on  plane  6  =  (J)*  or  $  =  0.  The  parameter  Sn 
has  the  value  of  +1  and  -1  on  the  upper  and  lower  fin  surface  respectively,  ks 
shown  in  Fig.  B-2,  the  radial  locations  of  the  inner  and  outer  fin  edges  are 
denoted  by  rg  and  r^,  and  the  quantities  (J)|  and  represent  the  angular 
locations  of  the  upper  and  lower  surfaces  of  the  fins  number  i  and  i+1 
respectively.  In  the  above,  represents  the  fin  thickness,  which  for  the 

purposes  of  force  and  moment  calculations  has  not  been  set  to  zero.  For  test 
cases  discussed  in  Section  5,  c()  is  set  to  4)^  and  hence  fin  thickness  is  neglected. 
In  the  case  of  pitch  plane  symmetry: 

III  =  =  3Mx  ^  Q 

3z  3z 

The  integrals  appear xrg  ii.  (B-1)  to  (B-6)  are  computed  numerically  at  each 
step,  z  =  Z^.  The  body  integrals  are  all  of  the  form: 

-F  <l>^ 


>?  =  2  /  V  ((t))d(j) 


These  are  evaluated  using 


ip  =  i' 


where:  ~  f  v((t))d4)  ; 


0  5  ,  1  =  <*>* 


In  evaluating  ip,  the  body  surface  pressure  on  fin  planes  is  taken  to  be  an 
average  of  the  surface  pressure  on  the  upper  and  lower  body  fin  junctions.  It 
is  more  convenient  to  perform  the  integration  of  equations  in  the  computational 
plane.  In  the  symmetric  case  the  above  becomes: 


li;  -  2  r  =  2/^  v(Y)  dY  . 

This  is  numerically  integrated  on  the  uniform  computational  mesh  using  Simpson's 
rule  in  the  form 

2AY 

^  +  4v2  +  2Vj  +  4v^  +  2V5  +  .  . 


B-2 
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and 


i)(  =  ^ ^  +  4v2  +  2V3  +  4v4  +  2V5  +  ... 

..  +  2vj^_3  +  4v^j_2  +  ^  (5V}4_3  +  3vj^)],  if  M  is  even. 


In  the  above,  =  v(Yn,)  ,  =  1  and  AY  =  l/(M-l)o.  In  the  last 

expression,  the  trapezoidal  rule  is  used  for  the  subinterval  [Yj^_3,  .  In 

the  nonsymmetric  problem  i'  is  written  as 


In  this  case,  the  Integrands  are  periodic  functions  of  Y  with  period  1  (l.e., 
v(0)  =  v(l))  and  the  Simpson's  rule  becomes 

ip  =  AX  [2^2  ^3  ^'*4 


.  .  +  2vj^_2  +  4vj|  ],  if  M  is  odd 


and 


'I' 


AX  [6v2  +  3v2  +  2V3  +  4v^  +  2V5  +  ... 


..  +  2vj^_2  +  ^''M-ll  >  M  is  even, 

where  =  v(Yni),  =  1-AY  and  AY  =  1/M. 

Integration  over  the  fin  surfaces  is  more  cumbersome  since  inner  and  outer 
fin  edges  in  every  crossflow  plane  need  not  coincide  with  a  grid  point.  Also 
the  fin  may  consist  of  many  or  only  a  single  grid  point.  In  order  to  allow  for 
leading  edge  functions  which  are  double  valued  in  z,  such  as  on  an  arrow  fin, 
the  fin  is  not  required  to  extend  to  the  body  surface.  Fin  surface  integrals 
have  the  form: 


f  -  J"  ^(r)dr 


The  integration  is  carried  out  in  computational  coordinates  and  leads  to: 


y 


%)dX 


X 

f  v(X)dX 
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Here  Xj  and  Xq  are  the  inner  and  outer  fin  edge  locations.  Consider  the 
case  where  the  inner-most  and  outer-most  grid  points  on  the  surface  of  the  fin 

occur  at  Xj,  and  respectively.  The  fin  surface  integral  is  expressed  in 
1  o 

three  parts: 


-V 

r"!  po  r° 

•¥  =  I  vdX  +  J  vdX  +  J  vdX 


(B-7) 


rhe  first  and  last  terms  are  approximated  by: 


/Li^  X 

VdX  =  v^  (X^  -  X.)  ;  J*°  VdX  =  v^^^  (Xo  -  X^  ) 


where  Vj^  =  v(Xri  )  >  '^n  °°  ''(^n  )•  remaining  term  is  evaluated  using 

1  *  i  o  o 

Simpson's  rule  which  requires  that  there  be  an  odd  number  of  points  on  the  fin. 
If  the  fin  contains  an  even  number  of  points,  the  outer  most  interval 
[Xq  -  AX,  Xq  ]  is  evaluated  using  the  trapezoidal  rule  while  the  remaining 
i  i 

interval,  which  now  consists  of  an  odd  number  of  points  is  evaluated  using 
Simpson's  rule.  This  leads  to  the  following  integration  formula  for  terms  1,  2, 
and  3 : 


/  VdX  =  V„  (Xq  -  Xq  -  +  Vq  (Xq  >  X^  + 


+1  +  2Vq  ^.2  4Vq_^.3 


4  I  M 

no-lj  3 


/  VdX  =  Vq  (^  +  .X^  -  Xq  )  +  Vq  _i  Vq_(Xq  -  .Xi  ^  + 

7  n  H  rt  < 


(^'’n.+l  +  2Vq  .j.2  +  4v  +3 
I  1  1  1 


...  4v 


n  _2  iM  -1  4) 

o  J  3  o  2 


where  AX  =  1/N 

The  code  also  computes,  at  each  computational  step,  the  force  and 

moment  vectors  acting  on  the  body  truncated  at  z  =  Z^.  These  quantities  are 
defined,  for  example,  by 


B-4 
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u  7  a  r 

Fa(Z  )  -  /  - S  dz  (B-8) 

0  ®  ^ 

with  similar  expressions  for  the  other  truncated  force  and  moment  components* 

The  Integrals  of  the  type  (B-8)  are  evaluated  numerically  using  the  trapezoidal 
rule;  i.e.. 


k+1  k 

Fa(Z  )  =  Fa(Z  )  + 


z=«Z 


k  + 


3F 

3Z 


z”Z 


k+1 


with  similar  expressions  for  the  other  force  and  moment  coefficients.  Note  that 
this  calculation  requires  the  force  and  moment  on  the  body  truncated  at  the 
initial  plane  z  =  Zq.  These  quantities  must  be  given  along  with  the  initial 
flow  field-  data. 


The  final  results  are  presented  in  coefficient  form  by  dividing  the  force 

P_  , 

components  and  their  derivatives  by  —  V*  A  and  the  moment  components  and 

2  <»  ref 

their  derivative  by  £*  V^(A  )  (z  ),  where  A  and  Z  are  reference  area  and 

2  <»  ref  ref  ref  ref 

reference  lengths  respectively. 

The  centers  of  pressure  in  the  pitch  and  yaw  planes  are  also  calculated 
using: 


for  Fy  f  0  and  F^  0. 


(Zcp)p  “  [Zc  + 

<Zcp)y  *  [Zc  “  V^y] 


1 

S 
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